Penalty Stiffness

The respective penalty stiffness *kP*, present in (66), is defined in the directions coincident with the constraint directions of which the first two are tangent and the third one is perpendicular to the mid-surface. The stiffness form corresponds to the displacement dofs that represent the mentioned sums and differences of the top and bottom dofs. The definition of the penalty stiffness reads:

*e*

$$
\stackrel{\varepsilon}{k}\_P = \int\_{R\_\varepsilon} \stackrel{\varepsilon}{\mathbf{N}}^T \mathfrak{n} P \mathfrak{n}^T \stackrel{\varepsilon}{\mathbf{N}} \coth(\mathfrak{f}) \, d\xi\_3 d\eta\_i. \tag{75}
$$

The penalty parameter has to be larger than the terms of *e kM* and the terms of *e kC*, i.e., *P* max{(*kM*)*ij*; *i*, *j* = 1, 2, ..., *nM*}, where (*kM*)*ij* represent terms of the stiffness matrix and *nM* stands for the corresponding number of displacement dofs. Additionally, one has *P*  max{(*kC*)*ij*; *i* = 1, 2, ..., *nM*; *j* = 1, 2, ..., *nE*}, with (*kC*)*ij* being terms of the piezoelectricity matrix and *nE* denoting the corresponding number of electric potential dofs.

The above form (75) of the penalty matrix is suitable for the geometry I of the classical and modified transition elements. In the above equation, *Re* stands for the element side being the subset of the boundary *R* between the first-order and transition model: *Re* ⊂ *R*. The normal vector *n* is defined as perpendicular to the mid-surface *Sm*, i.e., *n* ⊥ *Sm*. It also approximately satisfies the condition *n* ∈ *Re* – approximately due to the finite element approximation of the element geometry. Note that two tangential directions and the third normal direction are coincident with the local directions present in the constraint Equation (27). On the element level, they are defined by the relations (49) and (50). Note also that the presence of the vector *n* in the above definition results in imposition of the penalty terms in the third (normal) direction only, i.e., the constraints are applied in this direction only. The shape function matrix terms are non-zero on the boundary *Re*. Hence, the non-zero penalty terms appear for the element nodes placed on this boundary.

In the case of the classical and modified transition elements of the geometry II, the above definition has to take into consideration the fact that the constraints (27) are imposed on the element vertical edge *Le* only. The edge can be characterized with: *Le* ⊂ *R*. The normal vector is defined analogously as in the case of the geometry I, i.e., *n* ⊥ *Sm* and approximately *n* ∈ *Le*. With these remarks in mind, one can write

$$\stackrel{\mathcal{c}}{\dot{\mathbf{k}}}\_{P} = \int\_{L\_{\mathbf{r}}} \stackrel{\mathcal{c}}{\mathbf{N}}^{T} \mathfrak{n} P \mathfrak{n}^{T} \stackrel{\mathcal{c}}{\mathbf{N}} \mathrm{der}(f) \, d\xi\_{\mathcal{I}\mathcal{N}} \tag{76}$$

where, due to one dimensional integration, the value of der(*J*) is applied for transformation of the global directions into the normalized direction *ξ*3. As previously, the constraints concern the differences of the top and bottom displacement dofs. The non-zero penalty terms appear for the element nodes and the corresponding dofs placed on the boundary *Le*, in the transverse direction only.

#### 4.8.2. The Enhanced Transition Element

Two types of constraints are present in the case of the corresponding transition model. The first type are the gradually changing constraints of no elongation of the normals to the mid-surface, represented by the last but one term of (34). The considered type of the element requires also the additional constraints corresponding to the gradual change of the transverse order of approximation from *q* ≡ *I* = 1 to an arbitrary value *q* ≡ *I*, where *I* is the hierarchical model order, and represented by the last term of (34). As the constraints of both types are fully active on the boundary with the first-order model (*α* = 1) and gradually switched off (*α* → 0) towards the boundary with the hierarchical or three-dimensional model, the nodal (collocation) values of the complement of the switching functions to unity 1 − *α* will appear in the penalty stiffness definitions so as to activate the penalty terms.

### The Modified Stiffness Matrix

The relations (66)–(68) are still valid for the enhanced transition elements and constraints of no elongation of the normals to the mid-surface *Sm*. However, the degrees of freedom *k* and *l* corresponding to the top and bottom displacements, and the other dofs *m*, as well, are associated with different nodes now because the constraints of no elongation of the lines normal to the mid-surface are applied in the entire element volume *Ve*. In the case of the geometry I, the top dofs *k* are associated with the nodes *a*4, *a*5, *a*6, *a*10, *a*11, *a*12, *a*16 the bottom dofs *l* with the nodes *a*1, *a*2, *a*3, *a*7, *a*8, *a*9, *a*15, and the other dofs *m* with the nodes *a*13, *a*14, *a*17, *a*18, *a*19, *a*20. The geometry II needs the following assignments: the top dofs *k* correspond to the nodes *a*4, *a*5, *a*6, *a*10, *a*11, *a*12, *a*15, the bottom dofs *l* to the nodes *a*1, *a*2, *a*3, *a*7, *a*8, *a*9, *a*14, and the other dofs *m* to the nodes *a*13, *a*16, *a*17, *a*18. These assignments

hold also for the constraints related to the change of the transverse approximation order. In this case, the stiffness matrix modification is described with the relation

$$
\stackrel{\epsilon}{k}\_{M} = \stackrel{\epsilon}{k}\_{M} + \mathbf{Z}^{T} \stackrel{\epsilon}{k}\_{P} \mathbf{Z}. \tag{77}
$$

*e*

*e*

Comparing (66) and (77), one can notice that there are no matrices *R* and *R*−<sup>1</sup> present in the above relation. This means that the replacement of the top and bottom dofs with their sums and differences and inversely is not necessary for the constraints of the second type.

## The Modifying Operators

In the case of the constraints of no elongation of the normals, the definitions of the operators *e R* and *Z*, (69)–(71) and (72)–(79), respectively, are still valid. However, the new assignments of the top *k*, bottom *l*, and other *m* dofs to the element nodes for geometries I and II have to be taken into account while calculating the corresponding operator matrices. In the case of the constraints of the second type, which concern the other dofs *m* only, the following definition is appropriate

$$
\stackrel{\epsilon}{\mathbf{Z}} = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{Z}^2 \end{bmatrix}. \tag{78}
$$

As it can be seen, the first component block of (78) is *Z*<sup>1</sup> ≡ **0**. The second component block has the form

$$\mathbf{Z}^{\bar{2}} = \text{diag}\left[ \dots , \mathbf{Z}\_{\mathfrak{m}} , \dots \right], \quad \mathbf{Z}\_{\mathfrak{m}} = \text{diag}\left[ 1, 1, 1 \right], \tag{79}$$

which results in appearance of the penalty terms for all dofs *m* in all three global directions.

## The Penalty Matrices

The constraints of no elongation of the normals require the following definition of the penalty stiffness matrix:

$$\stackrel{\epsilon}{k}\_P = \int\_0^1 \int\_0^1 \int\_0^{-\frac{\pi}{2}2 + 1} \stackrel{\epsilon}{N}(I - \mathfrak{a}^T) \mathfrak{n} P \mathfrak{n}^T (I - \mathfrak{a}) \stackrel{\epsilon}{N} \det(f) \, d\xi\_1 d\xi\_2 d\xi\_3 \tag{80}$$

where normal vector is *n* ⊥ *Sm*, with *Sm* being the mid-surface within the symmetricthickness geometry. The presence of this vector results in application of the constraints in the normal direction only. In turn, the presence of the switching function *I*−*<sup>α</sup>* guarantees the gradual change of the penalty terms. It employs the identity matrix *I* and the matrix of nodal values of the blending function:

$$\mathfrak{a} = [\ldots \mathfrak{a}\_m \ldots], \quad \mathfrak{a}^T = [\ldots \mathfrak{a}\_m^T \ldots]^T,\tag{81}$$

and the nodal blocks *αm* are diagonal and correspond to the dofs *m*, other than the top and bottom ones:

$$\mathfrak{a}\_{\rm III} = \begin{bmatrix} \mathfrak{a}\_{\rm III} & 0 & 0 \\ 0 & \mathfrak{a}\_{\rm III} & 0 \\ 0 & 0 & \mathfrak{a}\_{\rm III} \end{bmatrix},\tag{82}$$

and where *αm* ∈ 0, <sup>1</sup>. Its value depends on the longitudinal location of the node, the dof *m* is associated with. In other words, *αm* = *<sup>α</sup>*(*ξ*1,*<sup>m</sup>*, *ξ*2,*m*), where *ξ*1,*m* and *ξ*2,*m* represent first two normalized coordinates for the dof *m*.

In the case of the constraints of the second type one has

$$
\stackrel{\mathcal{c}}{\mathbf{k}}\_P = \int\_0^1 \int\_0^1 \int\_0^{-\frac{\mathbf{r}}{\xi\_2} + 1} \stackrel{\mathcal{c}}{\mathbf{N}}^T (I - \mathfrak{a}^T) P (I - \mathfrak{a}) \stackrel{\mathcal{c}}{\mathbf{N}} \det(f) \, d\xi\_1 d\xi\_2 d\xi\_3. \tag{83}
$$

These constraints also change gradually and are applied in all three global directions for degrees of freedom *m* of the element. The definitions (81) and (82) of the switching function remain unchanged.

It should stressed that the above definitions (80) and (83) are valid for both the transition elements acting between the piezoelectric elements based on the first-order shell theory and the elements based on either the hierarchical shell theory or the three-dimensional theory. The difference is that, in the first case, the entire surface *Sm* (characterized with *ξ*3 = 1 2 ) represents the shell mid-surface, while, in the second case, it represents the reference surface only (defined with *ξ*3 = 1 2 ). The reference surface becomes the mid-surface only on the boundary with the first-order model *R* ∩ *Sm*.

## **5. Numerical Tests**

Two aspects of performance of the proposed transition piezoelectric elements will be addressed in this section. The first one is the ability to remove high stress gradients between the 3D-based transition and 3D-based basic elements in the chosen model problems. The mechanical models of the basic elements conform to the 3D-based first-order shell theory and the 3D-based hierarchical shell or three-dimensional theory. The electric field is characterized with the same dielectric model conforming to either three-dimensional or 3D-based hierarchical symmetric-thickness theory.

The second aspect will be comparison of the solution convergence for the same model problems as in the first aspect with three versions of the transition elements employed. The model problems will concern the bending-dominated plate and shell and membranedominated shell. These problems manifest different and numerically demanding stress and strain states. Explanation of the differences between these three cases of the mechanical fields can be found in Reference [53]. The electric field and electric displacement states will be typical for piezoelectricity.

It should be stressed that our model problems presented below serve the numerical assessment of the presented algorithms. Because of that the loads and charges are assumed so as the mechanical and electric parts of the potential electro-mechanical energy (or potential co-energy) are of the same order – the case very demanding from the numerical point of view as the sign of the electro-mechanical potential energy may change in the convergence studies where the error of the energy is displayed as a function of the number of dofs in the problem. Note that, in the two technologically important actuating and sensing modes of piezoelectric operation, one part of the co-energy dominates the other one. Because of this, such modes are less demanding numerically.

#### *5.1. Model Problems and Methodology*

#### 5.1.1. Model Structures and Their Geometry

Three piezoelectric domains corresponding to the bending-dominated plate, bendingdominated half-cylindrical shell and a half of the membrane-dominated cylindrical shell are presented in Figures 5–7. The plate dimensions are *l* × *l*, the length of both shells is equal to *l*, while the circumferential length of a half of the cylinder is *l* ≈ 2*π R* with *R* being the radius of the mid-surface of the shells. All the structures are of the same thickness *t*. The following values are set in our tests: *l* = 3.1415 × 10−<sup>2</sup> m, *R* = 1.0 × 10−<sup>2</sup> m, *t* = 0.03 × 10−<sup>2</sup> m.

#### 5.1.2. Material, Load, Charge and Boundary Conditions

All three structures are made of the same piezoelectric material of material constants typical for piezoceramics [57]. The Young modulus is *E* = 0.5 × 10<sup>11</sup> N/m2, the Poisson ratio equals *ν* = 0.294, the dielectricity constant under constant stress is *δ* = 0.1593 × 10−<sup>7</sup> F/m, while the piezoelectric constants under constant stress are equal to *d*13 = *d*23 = −0.15 × 10−<sup>9</sup> C/N, *d*33 = 0.3 × 10−<sup>9</sup> C/N, and *d*52 = *d*61 = 0.5 × 10−<sup>9</sup> C/N.

The plate and bending-dominated shell are loaded with the vertical traction of the value *p* = −0.4 × 10<sup>6</sup> N/m2, while the membrane-dominated shell with the internal pressure of the same value *p* acting in the outward direction. All the structures are electrically charged on their upper (or outward) surface where the surface charge density is *c* = 0.2 × 10−<sup>1</sup> C/m2. On the lower (or inward) surfaces the charge density is set to 0.

**Figure 5.** The domain of the bending-dominated piezoelectric plate.

**Figure 6.** The domain of the bending-dominated piezoelectric half-cylindrical shell.

**Figure 7.** A half of the membrane-dominated piezoelectric shell domain.

The mechanical boundary conditions are as follows. The plate is clamped along its four edges. The bending-dominated shell has two straight edges clamped and two curved edges free. In the case of the membrane-dominated shell, there is no rotation along the curved edges. The electrical boundary conditions assume grounding along: four edges of the plate, along two straight and two curved edges of the half-cylindrical shell, and two curved edges of the cylindrical shell.

### 5.1.3. Discretization and Methodology

Due to symmetry of the geometry, loading, electric charge, and boundary conditions, only quarters of the bending-dominated domains and an octant of the membranedominated domain are taken into account in or numerical tests and presented in the figures below. The numerical settings are as follows. We apply the meshes 4 × 4 × 2 of the prismatic elements in the displayed quarter or octant parts of the structures in our tests. The longitudinal order of approximation within the mechanical and electric field is the same and equal *p* = *π* = 1, 2, ... , 8. The transverse orders of approximation of both fields are equal to *q* = *ρ* = 2 in the hierarchical shell zone of the domains, *q* = *ρ* = 1 in the first-order shell zone, and the changing (from 2 to 1) order in the transition zone constituting a layer composed of couples of prismatic transition elements. The number of the first-order and hierarchical (three-dimensional) elements changes from the entirely hierarchical shell structure to the shell structure composed of the first-order shell elements only. The intermediate (or mixed) case includes a single layer of piezoelectric transition elements. The exemplary cases for the plate and shells are presented in Figures 8–10.

**Figure 8.** Mixed models of the bending-dominated plate.

**Figure 9.** Mixed models of the bending-dominated shell.

**Figure 10.** Mixed models of the membrane-dominated shell.

The presented meshes are composed of one layer of hierarchical elements (yellow and marked as MI in the figure), one layer of the transition elements (green and marked as RM/MI), and two layers of the first-order elements (blue and marked as RM in the figure).

Our parametric studies are limited to changes of the longitudinal approximation orders *p* ≡ *π*, the ratios of the numbers of elements of various types, and the types of the transition elements for three model problems. The mentioned ratios are equal to *r* = 1.0, 0.5, 0.0 and represent the number of the first-order elements divided by the sum of the transition and hierarchical elements along the characteristic longitudinal dimensions of the structures.

#### *5.2. Ability to Remove High Stress Gradients*

In order to present the ability of three (classical, modified, and enhanced) types of the transition elements to remove the high stress gradients between the transition elements and the neighboring basic (hierarchical and first-order) elements, exemplary (*p* ≡ *π* = 6) distributions of the global normal stress *σ*33 (marked as *szz*) for three model problems and three types of the transition elements are presented. These distributions correspond to the mixed models from Figures 8–10. The stress distributions for the plate problem, with the classical, modified and transition elements employed, are presented in Figures 11–13, respectively. Figures 14–16 display the analogous three stress distributions for the bending-dominated shell. The last figures, Figures 17–19, correspond to the membranedominated problem with three types of the transition elements applied. Note that, in the figures, the stress values are written in numerical form not in analytical form due to capabilities of the numerical code used.

In the case of the plate and membrane problems, high solution gradients are visible when the classical transition elements are employed. These gradients are reduced for the mixed models with the modified and enhanced elements applied. In the case of the bent shell, high gradients are not visible for the classical elements as the displayed global normal stress component is a combination of the transverse and longitudinal stresses.

**Figure 11.** Transverse stress—the bent plate with the classical transition elements.

**Figure 12.** Transverse stress—the bent plate with the modified transition elements.

**Figure 13.** Transverse stress—the bent plate with the enhanced transition elements.

**Figure 14.** Third normal stress—the bent shell with the classical transition elements.

**Figure 15.** Third normal stress—the bent shell with the modified transition elements.

**Figure 16.** Third normal stress—the bent shell with the enhanced transition elements.

**Figure 17.** Third normal stress—the membrane with the classical transition elements.

**Figure 18.** Third normal stress—the membrane with the classical transition elements.

**Figure 19.** Third normal stress—the membrane with the enhanced transition elements.

In order to quantify the presented graphical results concerning stress gradients, the Table 1 is presented. Therein, gradients (jumps/differences) of the global stress *σ*33 between the face centers of the neighboring elements are displayed for three model structures. Five models for each of three structures are included: the entirely hierarchical (MI) shell model (*r* = 0.0), three mixed (RM/MI) models (*r* = 0.5) with the classical, modified, and enhanced (respectively, TR1, TR2, TR3) transition models, and entirely first-order (RM) model (*r* = 1.0). In the case of the bending-dominated (bent) plate, the chosen elements 31 and 24 correspond to 3D model and TR1 or TR2 or TR3 model. Both elements 24 and 23 correspond to either TR1 or TR2 or TR3 model, while the elements 23 and 16 to either TR1 or TR2 or TR3 model and RM model. In the case of the bending-dominated (bent) shell, the corresponding couples of elements are: 25 and 28, 28 and 27, and 27 and 30. In the case of the membrane (membrane-dominated shell), the corresponding couples are: 6 and 13, 13 and 14, 14 and 21.

**Table 1.** Summary of stress gradients [N/m2] for three model problems.


\* The presented result values correspond to *p* ≡ *π* = 6.

In the case of the plate, the third normal stress *σ*33 corresponds to the transverse normal stress. In the entirely hierarchical (MI) shell model, the observed values of gradients are small and due to finite element discretization. In the entirely first-order (RM) model, the jumps are extremely small (correspond to the numerical zero), as the plane stress assumption holds in all elements. In the case of the mixed TR1 model, the jump between the transition and first-order elements 23 and 16 results from the three-dimensional stress state and plane stress state in the neighboring elements. Such jumps are diminished to numerical zero for the cases of the TR2 and TR3 models, as the plane stress assumption holds on both sides of the neighboring faces of the transition and first-order elements. In both shell cases, the presented jumps of *σ*33 result from both discretization and plane stress assumption. The exact analysis for the transverse direction requires taking jumps for the stress components *σ*11 and *σ*13 = *σ*31 into account. Such an analysis (not revealed here) leads to exactly the same observations as in the case of the plate. For the membrane case, where the contribution of *σ*33 to the transverse normal stress is higher than the influence of two other stress components, the stress jumps in *σ*33 between elements 14 and 21 are smaller in the case of the TR3 model than in the cases of the TR2 and TR1 models. In the bent shell case, the jumps in *σ*33 between elements 27 and 30 grow from TR1 through TR2 to TR3, but this leads to zero transverse normal stress values in both neighboring elements for the TR2 and TR3 models. Demonstration of this needs inclusion of the stress components *σ*11 and *σ*13 = *σ*31 again.

#### *5.3. Solution Convergence of the Problems*

In the convergence studies, the logarithm of the approximation error in potential energy is presented as a function of the logarithm of the number of dofs *N* of the problem. The corresponding curves are *p*-convergence ones as the number of dofs in the problems increase due to the change of the longitudinal order of approximation *p* within elements. In order to present the solution convergence for the three model problems with three types of the transition elements applied, three threesomes of pictures are displayed. The first three pictures (Figures 20–22) correspond to the bending-dominated plate problem, the second threesome (Figures 23–25) to the bending-dominated shell problem, while the last three drawings (Figures 26–28) to the membrane-dominated shell. In the figures, the first, second and third drawings correspond to the cases of the classical, modified, and enhanced transition elements applied.

The mentioned approximation error is defined as the difference of two energy measures of which *E* corresponds to the numerical solution under consideration and the reference solution *Er* plays a role of the exact solution, namely:

$$
\langle \text{energy error} = |E\_r - E|. \tag{84}$$

The energy measures represent the sum of absolute values of the mechanical and electrical parts of the electro-mechanical potential energy, i.e.,

$$E = |B(\mathfrak{u}, \mathfrak{u}) - \mathbb{C}(\mathfrak{u}, \mathfrak{\phi})| + |-b(\mathfrak{\phi}, \mathfrak{\phi}) - \mathbb{C}(\mathfrak{\phi}, \mathfrak{u})|.\tag{85}$$

Due to stationarity of the solutions *u* and *φ*, resulting from (8), the potential energy can be expressed with the components of the electro-mechanical co-energy (strain, electric field, and coupling energies) corresponding to the bilinear forms *B*, *b* and *C* from (9)–(11). The solutions for displacements *u* and electric potential *φ* are approximated with use of the global nodal vectors of the displacement dofs *qhpq* and potential dofs *ϕhπρ* obtained from solution of the problem (35).

**Figure 20.** Convergence curves—the plate problem, the classical transition elements employed.

**Figure 21.** Convergence curves—the plate problem, the modified transition elements applied.

**Figure 22.** Convergence curves—the plate problem, the enhanced transition elements applied.

**Figure 23.** Convergence curves—the bending-dominated shell with the classical transition elements.

**Figure 24.** Convergence curves—the bending-dominated shell with the modified transition elements.

**Figure 25.** Convergence curves—the bending-dominated shell with the enhanced transition elements.

**Figure 26.** Convergence curves—the membrane shell with the classical transition elements.

**Figure 27.** Convergence curves—the membrane shell with the modified transition elements.

**Figure 28.** Convergence curves—the membrane shell with the enhanced transition elements.

The proposed measure of the approximation error was applied as it guarantees monotonicity of the energy approximation error with changing value of *p* ≡ *π*, while the simple sum of both parts of the energy (equal to the electro-mechanical co-energy) may not. Note that such monotonicity is a must in the case of error-controlled adaptivity suggested in Reference [11,12].

In order to calculate the errors, the exact values *Er* of the energy measures are replaced with their best numerical approximations obtained from the so-called over-killed meshes, i.e., the meshes for which the longitudinal order of approximation is equal to *p* = *π* = 9.

In the case of two bending-dominated problems the convergence for the enhanced models is slightly higher than for the modified transition elements employed in the mixed models. Additionally, the latter elements produce slightly better convergence than the classical ones. In the case of the membrane-dominated problem, only the enhanced transition elements deliver convergence curves between the curves obtained for the basic models. The confirmation is presented in the Table 2, where the absolute values log(*Er*−*<sup>E</sup>*) of the error and relative error values (*Er*−*<sup>E</sup>*)/*Er* (both for *p* = 8, i.e., *E* ≡ *E*8) are included. In addition, the averaged convergence rates are given and equal to <sup>−</sup>log(*δ*8/*δ*7)/log(*N*8/*N*7), *δ*8 = *Er*−*E*8, *δ*7 = *Er*−*E*7, where *E*8, *N*8 and *E*7, *N*7 represent the energies and numbers of dofs for *p* = 8 and *p* = 7.

**Table 2.** Error and convergence summary for three model problems.


\* The presented result values correspond to *p* ≡ *π* = 8.
