Displacement assumptions

The starting point for defining the kinematic assumption within this type of the transition model is the definition of the displacement field of the 3D-based hierarchical shell models [13,14]. This field can be expressed with use of the local or global displacements, *u* or *u*, which are related by the transformation matrix *<sup>θ</sup>T*, i.e., *u* = *θTu*. The local displacement components are equal to:

$$\begin{split} u\_j^{\prime\prime} &= \theta\_{\bar{\imath}\bar{\jmath}} u\_i = \sum\_{l=0}^{l} f\_l(\mathbf{x}\_3^{\prime}) \theta\_{\bar{\imath}\bar{\jmath}} u\_i^{l} = \sum\_{l=0}^{l} f\_l(\mathbf{x}\_3^{\prime}) u\_j^{\prime l} \\ &= \ f\_0(\mathbf{x}\_3^{\prime}) u\_j^{\prime 0} + \sum\_{l=1}^{l-1} f\_l(\mathbf{x}\_3^{\prime}) u\_j^{\prime l} + f\_l(\mathbf{x}\_3^{\prime}) u\_j^{\prime l} \\ &= \left(\frac{1}{2} - \frac{\mathbf{x}\_3^{\prime}}{t}\right) u\_j^{0} + \left(\frac{1}{2} + \frac{\mathbf{x}\_3^{\prime}}{t}\right) u\_j^{\prime l} + \sum\_{l=1}^{l-1} f\_l(\mathbf{x}\_3^{\prime}) u\_j^{\prime l} \\ &= \ \frac{1}{2} (u\_j^{0} + u\_j^{\prime l}) + \frac{\mathbf{x}\_3^{\prime}}{t} (u\_j^{\prime l} - u\_j^{\prime l}) + \sum\_{l=1}^{l-1} f\_l(\mathbf{x}\_3^{\prime}) u\_{j\prime}^{\prime l} \quad \mathbf{x} \in V\_{\mathbf{z}} \end{split} \tag{33}$$

where *V* is the symmetric-thickness domain. In the first line above, *ulj*, *l* = 1, 2, ..., *I*, are the functions describing changes of displacements in the longitudinal directions. In the case of the symmetric-thickness geometry, these function are the same for each *l*. The functions *fl*(*x*3), *l* = 1, 2, ..., *I* stand for the polynomial functions in the direction normal to the midsurface. Along this direction the coordinate *x*3 ∈ (− *t*2 , *t*2 ) is employed. The polynomial functions in this direction span the respective polynomial space of order *I*. Hence, the quantity *I* represents the hierarchical model order. It is worth mentioning that with *I* → ∞ one reaches three-dimensional elasticity description of the above displacement field.

Following our propositions from Reference [14,16], in the second and third lines of (33), we define the first and the last functions, *f*0(*x*3) and *fI*(*x*3), as linear ones and the rest functions *fl*(*x*3), *l* = 1, 2, ..., *I* − 1 as polynomials of order *I*. After some arithmetic manipulations, one can reach the final form (line) of (33).

In the case of the enhanced transition piezoelectric model, its transverse displacement field is assumed to change from the state of no elongation of the lines perpendicular to the mid-surface to the state of lack of such an assumption. In the same manner, the order of the transverse displacement field is assumed to change from *I* = 1 to an arbitrary value *I*. So as to implement such changes, the blending function *α*, now playing the role of a gradually switching function, is applied, namely:

$$u\_j^{\prime} = \frac{1}{2} (u\_j^{\prime 0} + u\_j^{\prime I}) + \frac{\mathbf{x}\_3^{\prime}}{t} (u\_k^{\prime 0} - u\_k^{\prime I}) + a \frac{\mathbf{x}\_3^{\prime}}{t} (u\_3^{\prime 0} - u\_3^{\prime I}) + a \sum\_{l=1}^{I-1} f\_l(\mathbf{x}\_3^{\prime}) u\_{j\prime}^{\prime l}, \quad \mathbf{x} \in V,\tag{34}$$

where *V* is the symmetric-thickness domain or the domain with the mid-surface defined on *R* at least, while *j* = 1, 2, 3 and *k* = 1, 2. The later index *k* corresponds to the first two, longitudinal directions, while the index 3 is used for the transverse direction. It can be noticed that with *α* ≡ 1, (34) becomes (33) which characterizes the hierarchical or three-dimensional models, while, for *α* ≡ 0, (34) becomes (25) and (26) which correspond to the first-order model. In this case, *l* = 0 and *l* = 1 refer to the bottom (*l* ≡ *b*) and top (*l* ≡ *t*) surfaces.

#### **4. Algorithms of the Transition Piezoelectric Elements**

The starting point for our finite element considerations is the variational Formulations (8)–(11). Derivation of the corresponding finite element formulation needs discretization of the domain *V* and into finite element subdomains *Ve* and introduction of polynomial interpolation of the displacements *u* and electric potential *φ* (see Reference [11,12]) with use of the global vectors of the nodal, displacement *qhpq* and potential *<sup>ϕ</sup>hπρ*, degrees of freedom (dofs). The applied interpolations are suited for the *hpq*- and *hπρ*-adaptive approximations of the displacements and potential fields, with *h* being the common element size for both fields and independent longitudinal and transverse approximation orders *p*, *q* and *π*, *ρ* within the mechanical and electric fields, respectively. All these lead to

$$\begin{aligned} \mathbf{K}\_M \boldsymbol{q}^{\mathrm{hp}} - \mathbf{K}\_{\mathbb{C}} \boldsymbol{q}^{\mathrm{h\tau\rho}} &= \boldsymbol{F}\_V + \mathbf{F}\_{\mathbb{S}\prime} \\ \mathbf{K}\_{\mathbb{C}}^T \boldsymbol{q}^{\mathrm{hp}} + \mathbf{K}\_E \boldsymbol{q}^{\mathrm{h\tau\rho}} &= \boldsymbol{F}\_{Q\prime} \end{aligned} \tag{35}$$

where the components of the first equation correspond to: the first equation (9), the second Equation (9) and the first Equation (11), while the components of the second line result from: the second definition (11), the first definition (10) and the second definition (10). The global (of the entire domain) terms *K M*, *KC*, *KE* represent the global stiffness matrix, the global characteristic matrix of piezoelectricity, and the global characteristic matrix of dielectricity. The global vectors of nodal mass and surface forces are denoted as *FV* and *FS*, while the global vector of the surface charges as *<sup>F</sup>Q*.

#### *4.1. Normalized Geometry of the Transition Elements*

For three (classical, modified, enhanced) transition piezoelectric models two principal geometries of the prismatic elements has to be introduced in order to be able to join the basic elements corresponding to the three-dimensional or hierarchical models with the first-order model of piezoelectricity.

The example of the first transition geometry I, presented in Figure 1, includes one vertical edge (vertex nodes: *a*1, *a*4) corresponding to the 3D-based first-order shell theory [13,54] within the mechanical field. The other part of this field corresponds to the 3D-based hierarchical shell [13,14] or three-dimensional elasticity models [13,53]. Above, the following linear vertex (*<sup>a</sup>*2, *a*3, *a*5, *a*6) nodes and generalized: horizontal mid-edge (*<sup>a</sup>*7,*a*8, ..., *a*12), vertical mid-edge (*<sup>a</sup>*13, *a*14), mid-base (*<sup>a</sup>*15, *a*16), mid-side (*<sup>a</sup>*17,*a*18,*<sup>a</sup>*19) and middle (*a*20) nodes are applied within the geometry and displacement fields. For the element I, the electric field corresponds to the hierarchical symmetric-thickness or three-dimensional dielectricity models [11,16]. The linear vertex (*b*1, *b*2,. . . , *b*6) nodes, the generalized: horizontal mid-edge (*b*7, *b*8,. . . , *b*12), vertical mid-edge (*b*13, *b*14, *b*15), mid-base (*b*16, *b*17), mid-side (*b*18, *b*19, *b*20) and middle (*b*21) nodes characterize the approximated geometry and electric potential fields of the element.

**Figure 1.** An exemplary principal geometry I of the transition elements.

It is worth mentioning that the generalized (varying-order) nodes of both types may contain multiple degrees of freedom along the edges, within the bases and sides, and in the interior of the element. The related finite element approximations within the proposed element are all based on the general rules proposed by Demkowicz and others [1,2].

In the example of the second transition geometry II, shown in Figure 2, one deals with two vertical edges (nodes *a*1, *a*4 and *a*2, *a*5) and the side (nodes *a*7, *a*10) between these edges, all corresponding to the first-order mechanical model. The rest of the element (nodes *a*3, *a*6, *a*8, *a*9, *a*11, *a*12, *a*13, *a*14, *a*15, *a*16, *a*17, *a*18) conforms to the hierarchical shell model or the three-dimensional model of elasticity. These mechanical models are coupled with the symmetric-thickness hierarchical or three-dimensional models of dielectricity. The same nodes *b*1, *b*2, ..., *b*21, as for the geometry I, are employed in this case.

**Figure 2.** An exemplary principal geometry II of the transition elements.

#### *4.2. Element Characteristic Matrices and Vectors*

Replacement of the global quantities introduced in (24) by the corresponding sums of the element (or local) quantities leads to

$$\begin{aligned} \sum\_{\varepsilon} \left( \overset{\varepsilon}{\dot{\mathbf{k}}}\_{M} \overset{\varepsilon}{\boldsymbol{q}}^{hp\boldsymbol{q}} \right) &= \overset{\varepsilon}{\dot{\mathbf{k}}}\_{\mathcal{C}} \overset{\varepsilon}{\boldsymbol{q}}^{hp\boldsymbol{\pi}\boldsymbol{\rho}} \right) &= \sum\_{\varepsilon} \left( \overset{\varepsilon}{\boldsymbol{f}}\_{V} + \overset{\varepsilon}{\boldsymbol{f}}\_{S} \right), \\ -\sum\_{\varepsilon} \left( \overset{\varepsilon}{\boldsymbol{k}}\_{\mathcal{C}} \overset{\varepsilon}{\boldsymbol{q}}^{hp\boldsymbol{\eta}} \right) &= \overset{\varepsilon}{\boldsymbol{k}}\_{E} \overset{\varepsilon}{\boldsymbol{q}}^{hp\boldsymbol{\pi}\boldsymbol{\rho}} \end{aligned} \tag{36}$$

where *e* is the element number, while *qhpq* and *ϕhπρ* stand for the element, displacements and potential, dofs vectors.

*e*

The element stiffness matrix is defined in the standard way, namely:

*e*

$$\stackrel{\epsilon}{\mathbf{k}}\_M = \int\_0^1 \int\_0^1 \int\_0^{-\frac{\pi}{2} \mathbf{2} + 1} \stackrel{\epsilon}{\mathbf{B}}^T \mathbf{D} \stackrel{\epsilon}{\mathbf{B}} \det(\mathbf{J}) \, d\xi\_1 d\xi\_2 d\xi\_3 \,\tag{37}$$

*e*

with *e B* relating strains *ε* with the element displacement dofs *qhpq* and, hence, named the element strain-displacement matrix. The matrix *D* represents the matrix of elasticities for the isotropy case under consideration, while the term det(*J*) is the determinant of the Jacobian matrix transforming the global coordinates *x* into normalized ones *ξ* = [*ξ*1, *ξ*2, *ξ*3]*<sup>T</sup>*. Note that the applied normalized prismatic geometry of the element *e* can be characterized

with *Ve* = {(*ξ*1, *ξ*2, *ξ*3) ∈ [0, −*ξ*<sup>2</sup> + 1] × [0, 1] × [0, <sup>1</sup>]}. In the case of the thick- or thin-walled geometry, the first two normalized directions coincide with the longitudinal shell directions and the third normalized direction is the same as the transverse shell direction.

The definition of the element characteristic matrix of dielectricity reads

$$\stackrel{\epsilon}{k}\_E = \int\_0^1 \int\_0^1 \int\_0^{-\frac{\pi}{6}z+1} \stackrel{\epsilon}{b}^T \gamma \stackrel{\epsilon}{b} \det(f) \, d\xi\_1 d\xi\_2 d\xi\_3 \tag{38}$$

*e*

where *e b* allows expressing the electric field *E* with the element electric potential dofs *<sup>ϕ</sup>hπρ*. It is called the element field-potential matrix. The matrix *γ* is the constitutive constant matrix of orthotropic dielectricity under constant strain.

The element characteristic matrix of piezoelectricity coupling the mechanical and electric fields is defined as follows:

$$\stackrel{\text{c}}{\mathbf{k}}\_{\text{C}} = \int\_{0}^{1} \int\_{0}^{1} \int\_{0}^{-\frac{\pi}{6}\mathbf{2} + 1} \stackrel{\text{c}}{\mathbf{B}}^{T} \stackrel{\text{c}}{\mathbf{C}} \stackrel{\text{c}}{\mathbf{b}} \det(\mathbf{J}) \, d\xi\_{1} d\xi\_{2} d\xi\_{3} \,\tag{39}$$

with *C* being the constitutive matrix of orthotropic piezoelectricity under constant strain.

The element nodal forces vector due to the mass loading *f* is defined in the standard way, i.e.,

$$\stackrel{c}{\dot{f}}\_V = \int\_0^1 \int\_0^1 \int\_0^{-\frac{\pi}{5}\_2 + 1} \stackrel{c}{\mathbf{N}}^T f \det(\mathbf{J}) \, d\xi\_1 d\xi\_2 d\xi\_3,\tag{40}$$

with *N* standing for the standard shape functions matrix corresponding to the vectorial displacement field.

In addition the element nodal forces due to the surface loading *p* is defined in the standard way, for the bases and lateral sides, respectively:

$$\begin{aligned} \stackrel{\epsilon}{f}\_S &= \int\_0^1 \int\_0^{-\frac{\epsilon}{f\_2} + 1} \stackrel{\epsilon}{N}^T \not{p \, \text{cof}(f)} \, d\_{\mathbb{S}1}^x d\_{\mathbb{S}2}^x \\ \stackrel{\epsilon}{f}\_S &= \int\_0^1 \int\_0^1 \stackrel{\epsilon}{N}^T \not{p \, \text{cof}(f)} \, d\eta\_i d\_{\mathbb{S}3}^x. \end{aligned} \tag{41}$$

Above, *ηi*, *i* = 1, 3 is identical to the directions of the triangle edges, i.e., either *ξ*1 or *ξ*2 or 1 − *ξ*1 − *ξ*2, while the coefficient cof(*J*) transforms the surface element of the part *dSe* of the domain surface *S* from the global coordinates to the normalized ones.

Finally, the element nodal charges due to the charge density *c* are equal to

$$\begin{aligned} \stackrel{\varepsilon}{f}\_Q &= \int\_0^1 \int\_0^{-\frac{\varepsilon}{f\_2} + 1} \stackrel{\varepsilon}{n}^T \operatorname{c\ cost}(f) \, d\xi\_1 d\xi\_2\\ \stackrel{\varepsilon}{f}\_Q &= \int\_0^1 \int\_0^1 \stackrel{\varepsilon}{n}^T \operatorname{c\ cost}(f) \, d\eta\_i d\xi\_3 \end{aligned} \tag{42}$$

for the element bases and sides, respectively, with *en* denoting the standard shape function matrix suitable for the finite element approximation of the scalar electric potential field.

In the next subsection, we will present the changes that have to be implemented in the above finite element formulation depending on the applied transition model.

## *4.3. Strain-Displacement Matrix*

*e*

Three cases of the classical, modified, and enhanced transition elements corresponding to the analogous transition models are of interest in this section.

#### 4.3.1. The Transition Element Based on the Classical Model

The strain-displacement matrix for the classical transition model is assumed in the form corresponding to the three-dimensional elasticity and the hierarchical shell models, for the three-dimensional and thick- or thin-walled geometries, respectively. This means that the three-dimensional strain and stress state is applied for this classical model. This state is characterized by

$$
\varepsilon = \stackrel{c}{\mathbf{B}} q^{hpq} \tag{43}
$$

*e*

and reflects the relation between the local strains *ε* and global nodal displacement dofs *e* Thelocalstraindirectionscoincidentwithfirst

*<sup>q</sup>hpq*. are two tangent (the two or two longitudinal) and one normal (the third or transverse) directions within the structure.

The strain-displacement matrix *e B* will also be denoted as *B*<sup>1</sup> ≡ *eB* in our further considerations. The form of this matrix is derived from the standard approach applied to five-parameter thick-shell elements of the second order [63]. Here, we adopt it to the needs of the 3D-based hierarchical shell model. This form reads:

$$\stackrel{\varepsilon}{\mathbf{B}} = [\dots, \mathbf{B}\_{k}, \dots], \quad \stackrel{\varepsilon}{\mathbf{B}}^T = [\dots, \mathbf{B}\_{k}^T, \dots]^T. \tag{44}$$

The above blocks *Bk* ≡ *B*1*k* of the strain-displacement matrix correspond to the element displacement dofs at the element nodes and are composed of the following components:

$$B\_k = \begin{bmatrix} B\_{k,11} & B\_{k,12} & B\_{k,13} \\ B\_{k,21} & B\_{k,22} & B\_{k,23} \\ B\_{k,31} & B\_{k,32} & B\_{k,33} \\ B\_{k,41} & B\_{k,42} & B\_{k,43} \\ B\_{k,51} & B\_{k,52} & B\_{k,53} \\ B\_{k,61} & B\_{k,62} & B\_{k,63} \end{bmatrix}. \tag{45}$$

Any term *Bk*,*mn*, *m* = 1, 2, .., 6, *n* = 1, 2, 3 of each block *Bk*, where *m* corresponds to each of six components of the strain vector and *n* to each global direction of the displacement dof vector, can be determined as equal to:

$$B\_{k,1n} = \beta\_{k,1}\theta\_{1n}; \quad B\_{k,4n} = \beta\_{k,1}\theta\_{2n} + \beta\_{k,2}\theta\_{1n}$$

$$B\_{k,2n} = \beta\_{k,2}\theta\_{2n}; \quad B\_{k,5n} = \beta\_{k,2}\theta\_{3n} + \beta\_{k,3}\theta\_{2n} \tag{46}$$

$$B\_{k,3n} = \beta\_{k,3}\theta\_{3n}; \quad B\_{k,6n} = \beta\_{k,3}\theta\_{1n} + \beta\_{k,1}\theta\_{3n}$$

with the terms *<sup>θ</sup>pn* representing components of the transformation matrix *θ* which converts the local Cartesian directions *p* = 1, 2, 3 of the displacement dofs into the global ones. The terms *β<sup>k</sup>*,*<sup>p</sup>* can be calculated from

$$
\beta\_{k,p} = P\_{k,1}A^{1p} + P\_{k,2}A^{2p} + P\_{k,3}A^{3p}, \tag{47}
$$

where the quantities *Pk*,*<sup>r</sup>* = *∂Nk*/*∂ξr* are the partial derivatives, with respect to the normalized directions *r* = 1, 2, 3, of the shape function *Nk* from the corresponding shape function matrix. The auxiliary quantities *Arp* are equal to:

$$A^{rp} = j^{r1}\theta\_{1p} + j^{r2}\theta\_{2p} + j^{r3}\theta\_{3p\prime} \tag{48}$$

with the terms *jrn* representing terms of the inverse Jacobian matrix *J*−<sup>1</sup> and *<sup>θ</sup>np* being the components of the transposed transformation matrix *θT*. The latter matrix transforms derivatives with respect to the global Cartesian directions into derivatives with respect to the local Cartesian directions. The terms of the transformation matrices can be obtained from:

$$
\boldsymbol{\theta} = \begin{bmatrix} \theta\_{11} & \theta\_{12} & \theta\_{13} \\ \theta\_{21} & \theta\_{22} & \theta\_{23} \\ \theta\_{31} & \theta\_{32} & \theta\_{33} \end{bmatrix} = \begin{bmatrix} u\_1 & v\_1 & w\_1 \\ u\_2 & v\_2 & w\_2 \\ u\_3 & v\_3 & w\_3 \end{bmatrix} \tag{49}
$$

where the unit vectors components *up* = *Up*/*U*, *vp* = *Vp*/*V*, *wp* = *Wp*/*W*, *p* = 1, 2, 3 can be obtained from the vectors *U*, *V* and *W* equal to:

$$\mathbf{W} = \begin{bmatrix} \left\| \mathbf{x}\_1 / \partial \mathbf{\xi}\_1^x \\\\ \left\| \mathbf{x}\_2 / \partial \mathbf{\xi}\_1^x \right\| \\\\ \left\| \mathbf{x}\_3 / \partial \mathbf{\xi}\_1^x \right\| \end{bmatrix} \times \begin{bmatrix} \left\| \mathbf{x}\_1 / \partial \mathbf{\xi}\_2^x \\\\ \left\| \mathbf{x}\_2 / \partial \mathbf{\xi}\_2^x \right\| \\\\ \left\| \mathbf{x}\_3 / \partial \mathbf{\xi}\_2^x \right\| \end{bmatrix}, \quad \mathbf{U} = \begin{bmatrix} \left\| \mathbf{x}\_1 / \partial \mathbf{\xi}\_1^x \\\\ \left\| \mathbf{x}\_2 / \partial \mathbf{\xi}\_1^x \right\| \\\\ \left\| \mathbf{x}\_3 / \partial \mathbf{\xi}\_1^x \right\| \end{bmatrix}, \quad \mathbf{V} = \mathbf{W} \times \mathbf{U}, \tag{50}$$

where the partial derivatives *∂xn*/*∂ξr* = *jnr*, *n*,*r* = 1, 2, 3 are the terms of the Jacobian matrix *J*.

#### 4.3.2. The Modified and Enhanced Transition Elements

*e*

*e*

In the case of these two types of the transition element, the constitutive relations (29) and (30) hold within the coupled mechanical field. The strain and stress state changes from the three-dimensional one to the plane stress state. This can be expressed with the following relation

$$\varepsilon = \stackrel{\varepsilon}{B} q^{hpq} = \left[ a \stackrel{\varepsilon}{B}^1 + (1 - a) \stackrel{\varepsilon}{B}^2 \right] \stackrel{\varepsilon}{q}^{hpq} \, \tag{51}$$

where the components *B*<sup>1</sup> and *B*<sup>2</sup> correspond to the three-dimensional and plane stress states. The first component is defined in accordance with (44)–(50). The second component is divided into blocks corresponding to the displacement dofs at the element nodes, namely

$$\stackrel{c}{\mathbf{B}}^2 = [\dots, \mathbf{B}\_{k'}^2 \dots]. \tag{52}$$

Due to the plane stress assumption (21), any block *k* of the second component needs the following modification:

$$\mathbf{B}\_{k}^{2} = \begin{bmatrix} B\_{k,11} & B\_{k,12} & B\_{k,13} \\ B\_{k,21} & B\_{k,22} & B\_{k,23} \\ \frac{-\nu}{1-\nu} \left(B\_{k,11} + B\_{k,21}\right) & \frac{-\nu}{1-\nu} \left(B\_{k,12} + B\_{k,22}\right) & \frac{-\nu}{1-\nu} \left(B\_{k,13} + B\_{k,23}\right) \\ B\_{k,41} & B\_{k,42} & B\_{k,43} \\ B\_{k,51} & B\_{k,52} & B\_{k,53} \\ B\_{k,61} & B\_{k,62} & B\_{k,63} \end{bmatrix},\tag{53}$$

where the terms *Bk*,*mn*, *m* = 1, 2, 4, 5, 6, *n* = 1, 2, 3 are defined in accordance with (46)–(50). Because of the modified definition (51) of the strain vector *ε*, also any block *k* of the strain-displacementmatrixhastobemodifiedinaccordancewith

$$\mathbf{B}\_{k} = a\mathbf{B}\_{k}^{1} + (1 - a)\mathbf{B}\_{k}^{2} =$$

$$\begin{bmatrix} B\_{k,11} & B\_{k,12} & B\_{k,13} \\ B\_{k,21} & B\_{k,22} & B\_{k,23} \\ aB\_{k,31} + (1 - a)\frac{-\nu}{1 - \nu}(B\_{k,11} + B\_{k,21}) & aB\_{k,32} + (1 - a)\frac{-\nu}{1 - \nu}(B\_{k,12} + B\_{k,22}) & aB\_{k,33} + (1 - a)\frac{-\nu}{1 - \nu}(B\_{k,13} + B\_{k,23}) \\ B\_{k,41} & B\_{k,42} & B\_{k,43} \\ B\_{k,51} & B\_{k,52} & B\_{k,53} \\ B\_{k,61} & B\_{k,62} & B\_{k,63} \end{bmatrix} . \tag{54}$$

The blending function is defined in accordance with *α* = *<sup>α</sup>*(*ξ*1, *ξ*2, 12 ), i.e., as a function of the two normalized directions, *ξ*1 and *ξ*2, tangent to the mid-surface defined with *ξ*3 = 12 . Note that, for *α* ≡ 1 and *α* ≡ 0, the definition (54) transforms into (45) and (53), suitable for the basic (three-dimensional and first-order) piezoelectric elements, respectively.

## *4.4. Field-Potential Matrix*

For any type of the transition models: (16), (23), (24), and (30) plus (32), the definition of the electric field vector *E* remains unchanged. Because of this also the field-potential matrix has the same form for the classical, modified, and enhanced transition elements. This form corresponds to the three-dimensional or hierarchical symmetric-thickness dielectric models. For these models, the following relation holds:

$$E = \stackrel{\varepsilon}{b} \stackrel{\varepsilon}{\!\!\!\!\!\!\!\!\!\!\/} \,\tag{55}$$

*e*

Above, the locally defined electric field *E* is expressed with the scalar potential dofs *<sup>ϕ</sup>hπρ*. The form of the electric field-potential matrix under consideration is

$$\stackrel{\circ}{\mathbf{b}}^{\varepsilon} = [\dots, \mathbf{b}\_{l}, \dots]\_{\prime} \qquad \stackrel{\circ}{\mathbf{b}}^{T} = [\dots, \mathbf{b}\_{l}^{T}, \dots]^{T}\_{\prime} \tag{56}$$

where the blocks *bl* of the field-potential matrix correspond to the electric potential dofs of the element. The number of these blocks may be different to the number of the displacement dofs as the electric potential field is *hπρ*-interpolated, while the displacement field is *hpq*interpolated. The form of these blocks is:

$$\mathbf{b}\_{l} = \begin{bmatrix} b\_{l,1} \\ b\_{l,2} \\ b\_{l,3} \end{bmatrix}. \tag{57}$$

The terms *bl*,*p*, where *p* = 1, 2, 3 are the local Cartesian directions, are equal to

$$b\_{l,p} = p\_{l,1}A^{1p} + p\_{l,2}A^{2p} + p\_{l,3}A^{3p}.\tag{58}$$

Note that the quantities *pl*,*<sup>r</sup>* = *∂nl*/*∂ξ<sup>r</sup>*, *r* = 1, 2, 3 are the derivatives of the terms *nl* of the shape function matrix with respect to the normalized directions *r* = 1, 2, 3. The shape function *nl* corresponds to the electric potential dof *l* at the element node. The auxiliary terms *Arp*, *r*, *p* = 1, 2, 3 can be calculated from

$$A^{rp} = j^{r1}\theta\_{1p} + j^{r2}\theta\_{2p} + j^{r3}\theta\_{3p\prime} \tag{59}$$

i.e., analogously as in the case of the strain-displacement matrix. As before, the terms *<sup>θ</sup>np* belong to the transformation matrix *θ* which transforms the local Cartesian directions *p* = 1, 2, 3 into global Cartesian ones *n* = 1, 2, 3. The quantities *jrn* are the terms of the inverse Jacobian matrix *J*−<sup>1</sup> again.

### *4.5. Constitutive Constant Matrices*

In the case of the classical, modified, and enhanced transition elements, the constitutive, elasticity constant matrix from (37) takes the same form (13).

The classical transition element requires the dielectricity and piezoelectricity constant matrices, in (38) and (39), in the form defined by (18) and (19), respectively. As far as the modified and enhanced transition elements are concerned, one has to apply, in relations (38) and (39), the definitions included in (32) and (30), respectively.

## *4.6. Blending Functions*

As said above, blending functions describe the linear change of the stress state within transition elements from the three-dimensional to plane stress state. This change manifests itself by the presence of function *α* = *<sup>α</sup>*(*ξ*1, *ξ*2) in the strain definition of the modified and enhanced transition elements and in the piezoelectric and dielectric constant matrices for these two transition elements. Note that *ξ<sup>i</sup>*, *i* = 1, 2, 3 are the normalized directions within an element and *ξ*3 defines location of the reference surface within the element.

The blending function corresponding to the transition element of the geometry I is presented in Figure 3. The function is represented by the triangle spanned at the vertices *a*1, *a*5, *a*6 and takes the value of 0 for the vertical edge 1 corresponding to the mechanical field of the first-order in the transverse direction *ξ*3 and to the electric field of an arbitrary order in this direction. The blending function value of 1 corresponds to the mechanical and electric fields of arbitrary orders for the vertical edges 2 and 3. This function can be expressed with the affine coordinates *λ*1, *i* = 1, 2, 3 of the triangular element base or the normalized coordinates *ξ*1 and *ξ*2 as

$$a = 1 - \lambda\_i, \quad i = 1, 2, 3,\tag{60}$$

where *λ*1 = 1 − *ξ*1 − *ξ*2, *λ*2 = *ξ*1 and *λ*3 = *ξ*2.

**Figure 3.** The blending function for an example of the principal transition geometry I.

The geometry II of the transition element and the respective blending function are presented in Figure 4. The function is represented by the triangle spanned at the vertices *a*1, *a*2, *a*6 and equals 0 for two edges 1 and 2 corresponding to the mechanical field of the first-order in the transverse direction. This value is equal to 1 for the vertical edge 3 where the mechanical field is of arbitrary order in the transverse direction *ξ*3. The order of the electric field is arbitrary for all three vertical edges. The function can be expressed in the following way:

$$
\mathfrak{a} = \lambda\_{i\prime} \quad i = 1,2,3. \tag{61}
$$

**Figure 4.** The blending function for an example of the principal transition geometry II.

It should be noted that due to introduction of the linear blending functions into the vectors and matrices from (29)–(32), the polynomial order the integrands of (37)–(39) increases. This needs the respective increase of the number of the Gauss points in the numerical integration of these integrands.

### *4.7. Shape Function Matrices*

The shape function matrices *e N* corresponding to the displacement field of the transition piezoelectric elements are of the standard form for any type of the element. This form reads

$$\stackrel{\epsilon}{\mathbf{N}} = [\dots, \mathbf{N}\_{k'} \dots]\_{\prime} \quad \stackrel{\epsilon}{\mathbf{N}}^T = [\dots, \mathbf{N}\_{k'}^T \dots]^T \tag{62}$$

where the blocks *Nk* are diagonal and correspond to the vectorial displacement dofs *k* at the element node:

$$\mathbf{N}\_{k} = \begin{bmatrix} N\_{k} & 0 & 0 \\ 0 & N\_{k} & 0 \\ 0 & 0 & N\_{k} \end{bmatrix}. \tag{63}$$

In the case of the electric potential, the following standard form of the shape function matrix *en* is valid for the piezoelectric transition element of any type:

$$\stackrel{c}{\mathfrak{n}} = [\dots, \mathfrak{n}\_{l}, \dots], \quad \stackrel{c}{\mathfrak{n}}^T = [\dots, \mathfrak{n}\_{l}^T, \dots]^T,\tag{64}$$

where the blocks *nl*, corresponding to the scalar electric potential dof *l* at the element node, reduce to:

$$
\mathfrak{u}\_l = \begin{bmatrix} \ \mathfrak{u}\_l \end{bmatrix}.\tag{65}$$

It is worth noticing that the forms (44)–(45) and (64)–(65) corresponding to three proposed transition models are exactly the same as in the case of the three-dimensional and hierarchical models of piezoelectricity.

#### *4.8. Kinematic Constraints within the Elements*

In this section, the algorithms for imposition of the constraints of no elongation of the normals to the mid-surface, and for application of the constraints imposed on the generalized (varying-order) displacement dofs, as well, will be presented.

### 4.8.1. The Classical and Modified Elements Stiffness Matrix Modification

In this case the constraints of no elongation of the normals are assumed on the boundary *Re* between the first-order and transition elements. These constraints are imposed with the penalty method where the parameter *P* is employed. It represents the inverse of the penalty parameter and is a number sufficiently large to overwrite the stiffness matrix terms with the constraint terms and to replace the finite element equation with the constraint equation. Our approach takes advantage of the presentation of the method given in Reference [1,2]. Imposition of the constraints under consideration needs the following modification of the stiffness matrix:

$$
\stackrel{\varepsilon}{\dot{\mathbf{k}}}\_M = (\stackrel{\varepsilon}{\mathbf{R}}^{-1})^T \left( \stackrel{\varepsilon}{\mathbf{R}}^T \stackrel{\varepsilon}{\dot{\mathbf{k}}}\_M \stackrel{\varepsilon}{\mathbf{R}} + \stackrel{\varepsilon}{\mathbf{Z}}^T \stackrel{\varepsilon}{\dot{\mathbf{k}}}\_P \mathbf{Z} \right) \stackrel{\varepsilon}{\mathbf{R}}^{-1}, \tag{66}
$$

*e*

*e*

where *kP* is called the penalty matrix. Its terms include the mentioned parameter *P*.

Note that no modification of the force vectors *f V* + *f S* due to penalty method is required as one deals with zero constraints in the paper, e.g., the constraints of no elongation of the normals to the mid-surface. The same refers to other constraints presented in this paper.

*e*

The arithmetic operator matrices *e R* and *R*−<sup>1</sup> of the element *e* from (66) are capable of converting the top and bottom displacement dofs into their mean sums and differences, and conversely, respectively. The form of the first operator is related to the division of the element dofs into two groups: the pairs of top and bottom displacement dofs *q*1 and all other displacement dofs *q*2, namely:

$$q^{\mathcal{A}}{}^{\mu\eta} = \left[\begin{array}{c} q^1\\ q^2 \end{array}\right] . \tag{67}$$

where the two blocks consist of the following component sub-blocks:

$$\boldsymbol{q}^1 = [\ldots, \boldsymbol{q}\_{k'} \ldots, \boldsymbol{q}\_{1'} \ldots]^T, \quad \boldsymbol{q}^2 = [\ldots, \boldsymbol{q}\_{m'} \ldots]^T. \tag{68}$$

The top and bottom dofs, *k* and *l*, are associated with the nodes belonging to *Re* = *<sup>R</sup>*|*Ve* only. In the case of the transition element geometry I (Figure 3), the top dofs *k* are associated with the node *a*4, the bottom dofs *l* with the node *a*1, the other dofs *m* with the nodes *a*2, *a*3, *a*5, *a*6, *a*7, ..., *a*18. In the case of the geometry II, one has the following assignments: the top dofs *k* are associated with the nodes *a*4, *a*5, *a*10, the bottom dofs *l* with the nodes *a*1, *a*2, *a*7, the other dofs *m* with the nodes *a*3, *a*6, *a*8, *a*9, *a*11, *a*12, *a*13, ..., *a*20.
