1.1.3. Conclusions

In this paper, we will use piezoelectric transition elements combining the ideas present in the following works concerning the uncoupled problem of elasticity. Firstly, we will apply the 3D-based approach [13], where the solid [53], hierarchical shell [14], first-order shell [54], and transition [51] models are equipped with the three-dimensional dofs only. This approach will also be extended onto the dielectricity models [16]. We will also apply the idea of Dávila [35], where one deals with larger 3D part and smaller thin-walled part of the transition element. The stress state within the element will be either three-dimensional, as in Reference [35], or transition (from three-dimensional to plane one) [24]. The *hp*adaptive capabilities of the elements will be based solely on the ideas of Demkowicz and others [2]. Our 3D-based hierarchical finite element models will resemble the ideas presented, respectively, in Reference [55,56] for the cases of the three-dimensional and conventional (mid-surface) dofs.

#### *1.2. The Scope and Novelty of The Paper*

This paper concerns the algorithm of the solid-to-shell and shell-to-shell transition piezoelectric elements for adaptive analysis of electro-mechanical systems. Such systems can be composed of piezoelectric transducers and elastic structural elements, both of arbitrary shapes. In addition, effectivity of application of the proposed elements in such an analysis is of our interest in the paper. The elements are assigned for joining basic elements, namely piezoelectric shell elements of the first-order mechanical field and arbitrary-order electric field and piezoelectric elements of the hierarchical shell model (second- or higherorder) within the mechanical field and arbitrary-order within the electric field. Alternatively, the latter model can be replaced with the piezoelectric description based on hierarchical models of three-dimensional elasticity and three-dimensional dielectricity.

In this work, three variants of the transition models are applied. The first one, called classical, assures continuity of displacements between the basic models and continuity of electric potential between these models, as well. The second, modified transition piezoelectric model guarantees additional continuity of the stress field between the basic models. The third, enhanced model additionally enables continuous change of the strain state between the basic models. All three models are 3D-based ones, namely only threedimensional degrees of freedom are applied within the mechanical and electric fields. Based on the mentioned models, three types of the corresponding transition finite elements are introduced. The applied approximations are *hpq*/*hp*-adaptive ones, with *h* being the element size parameter, *p* standing for the longitudinal order of approximation, and *q* denoting the transverse order of approximation.

Due to basic character of the presented research, numerical effectiveness of these models and their finite element approximations is limited to two crucial aspects, namely ability to remove high stress gradients between the basic and transition models, and convergence

of the numerical solutions for the model problems of piezoelectrics. These solutions are obtained with and without the proposed transition elements. Error estimation and adaptivity control matters are out of scope of this work and will be presented in a separate paper.

Novelty of the paper results from the above unique features of the proposed transition elements and the orignal results concerning their effectiveness.

#### **2. Model Problems of Stationary Piezoelectricity**

First, we present the strong (or local) formulation of the model problem of stationary piezoelectricity. Due to further application of the matrix notation, convenient for finite elements, this formulation is expressed in the matrix form.

The mechanical equilibrium equations and geometrical (kinematic) relations read

$$\operatorname{div}\underline{\sigma} + f = \mathbf{0}$$

$$\underline{\mathfrak{e}} = \frac{1}{2} [\mathfrak{u} \otimes \operatorname{grad} + (\mathfrak{u} \otimes \operatorname{grad})^T] \quad \text{ $\mathfrak{x} \in V$ }.\tag{1}$$

where the matrix representations *σ* and *ε* of the symmetric, stress, and strain tensors, *σ*˜ and *ε* ˜, of the second rank are employed. Both relations are valid for each point *x* of the volume domain *V* of the piezoelectric. The vector *f* represents given mass loading vector, while *u* is the unknown vector of (mechanical) displacements.

Within the electric field, the electric equilibrium equation, namely the Gauss law with the assumed zero volume charges, holds together with the electric field vector *E* definition, i.e.,

$$\begin{aligned} \operatorname{div} d &= 0 \\ E &= -\operatorname{grad} \phi \end{aligned} \Big|\_{\prime} \quad \mathbf{x} \in V. \tag{2}$$

Above, the vector of electric displacements *d*, resulting from electrical charging of the piezoelectric, and the unknown scalar electric potential *φ* appear.

In order to express the stress tensor by the strain tensor, and the electric displacement vector by the electric potential one, the constitutive relations have to be added to the systems (1) and (2). The general form of these relations is consistent with the work of Reference [57]. The matrix form of them corresponds to the isotropic elasticity, orthotropic dielectricity, and orthotropic piezoelectricity. Because of this, the general, constitutive tensors of the elasticity, dielectricity, and piezoelectricity, *D* ˜ , *γ* ˜ , *C* ˜ , of the fourth, second, and third ranks, respectively, can be represented by the matrices *D*, *γ*, and *C*, of the sizes: 6 × 6, 3 × 3, and 6 × 3, respectively, i.e.,

$$\begin{cases} \sigma = \mathsf{D}\varepsilon - \mathsf{C}\mathsf{E} \\ d = \mathsf{C}^{T}\varepsilon + \gamma \mathsf{E} \end{cases} \right\}, \quad \mathsf{x} \in V\_{\prime} \tag{3}$$

where the coupling of the mechanical and electric fields results from the presence of the matrix *C*, while *σ* and *ε* stand for the six-component vectorial representations of the stress and strain tensors *σ*˜ and *ε*˜.

The above set of the partial differential, governing equations of piezoelectricity (1)–(3) has to be completed with the boundary conditions. The Neumann and Dirichlet boundary conditions, concerning, respectively, unknown functions derivatives and unknown functions themselves, become stress and displacement boundary conditions in the mechanical case, namely:

$$
\underline{\sigma}\mathfrak{n} = \mathfrak{p}, \quad \mathfrak{x} \in P \tag{4}
$$

$$
\mathfrak{u} = \mathfrak{w}\_{\prime} \quad \mathfrak{x} \in \mathbb{W}\_{\prime} \tag{5}
$$

where *p* and *w* are the given stress and displacement vectors, respectively, while *P* and *W* denote the parts of the surface *S* ≡ *∂V* of the piezoelectric body *V*, with the given values of the mentioned vectors. The vector *n* represents the unit vector, outward and normal to the boundary.

In the case of the electric field, the Neumann boundary conditions express equality of the internal charges (defined with the unknown electric displacements *d*) and given external charges of the surface density *c*, while the Dirichlet boundary conditions describe the given values *ω* of the unknown field of electric potential *φ*, i.e.,

$$d^I \mathfrak{n} = -\mathfrak{c}, \quad \mathfrak{x} \in Q,\tag{6}$$

$$
\phi = \omega, \quad \mathfrak{x} \in F. \tag{7}
$$

Above, *Q* and *F* represent the parts of the surface *S* of the piezoelectric domain *V*, where the electric charge and electric potential are known, respectively.

It should be mentioned that the general constitutive relations (3) allow any piezoelectric model based on coupling of the linear elastic and linear dielectric models. In our proposition, we apply the 3D-based hierarchical piezoelectric models proposed in Reference [15,58] and elucidated in Reference [11,12]. These models employ three-dimensional degrees of freedom within the electric and mechanical fields. The 3D-based elastic models may represent: three-dimensional elasticity, hierarchical and first-order shell models, as well as the solid-to-shell and shell-to-shell transition models. The 3D-based dielectric models include: three-dimensional dielectricity model and the hierarchical symmetricthickness models. These two groups of models were introduced and thoroughly described in Reference [13,14] and Reference [16], respectively. Note that the applied 3D-based approach allows analysis of complex piezoelectric structures within the above Formulations (1)–(3). In such structures, more than one piezoelectric model meet, including transition ones.

Thanks to the reciprocity theorem [59] of linear piezoelectricity, the above strong Formulations (1)–(7) can be converted into the corresponding variational (or weak) formulation. For the model problems of linear stationary piezoelectricity considered in this paper, the weak formulation reads

$$\begin{cases} \begin{array}{rcl} B(\boldsymbol{\upsilon},\boldsymbol{\mu}) - \mathbb{C}(\boldsymbol{\upsilon},\boldsymbol{\phi}) &=& L(\boldsymbol{\upsilon}) \\ -\mathbb{C}(\boldsymbol{\psi},\boldsymbol{\mu}) - b(\boldsymbol{\psi},\boldsymbol{\phi}) &=& -l(\boldsymbol{\psi}) \end{array} \end{cases} \tag{8}$$

for all admissible displacements *v* ∈ V and all admissible potentials *ψ* ∈ Ψ, with the space definitions: V = {*v* ∈ (*H*<sup>1</sup>(*V*))<sup>3</sup> : *v* = **0** on *W*} and Ψ = {*ψ* ∈ *H*<sup>1</sup>(*V*) : *ψ* = 0 on *<sup>F</sup>*}. The solution functions belong to the spaces: *u* ∈ *w* + V and *φ* ∈ *ω* + Ψ, where the given values of *w* and *ω*, consistent with the Dirichlet boundary conditions (5) and (7), are called lifts (compare Reference [1]) within the displacements field and potential field, respectively.

The bilinear and linear forms of the first line of the relation (8) denote the virtual strain energy of the elastic field characterized with the matrix *D* of elasticities and the virtual work of the external body *f* and surface *p* loadings. The following definitions of these forms hold: 

$$\begin{aligned} \mathbf{B}(\boldsymbol{\mathfrak{v}}, \boldsymbol{\mathfrak{u}}) &= \int\_{V} \boldsymbol{\mathfrak{e}}^{T}(\boldsymbol{\mathfrak{v}}) \mathbf{D} \, \boldsymbol{\mathfrak{e}}(\boldsymbol{\mathfrak{u}}) \, dV \\ \boldsymbol{L}(\boldsymbol{\mathfrak{v}}) &= \int\_{V} \boldsymbol{\mathfrak{v}}^{T} \boldsymbol{f} \, dV + \int\_{P} \boldsymbol{\mathfrak{v}}^{T} \boldsymbol{\mathfrak{p}} \, dS \end{aligned} \tag{9}$$

In the second Equation (8), the bilinear form represents the virtual energy of the electric field characterized with the dielectricity constant matrix *δ* under constant strain, while the linear form is equal to the virtual work of the surface charges of the density *c*, namely

$$\begin{aligned} \label{eq:SDAR} \dot{\psi}(\psi, \phi) &= \int\_V E^T(\psi) \, \delta E(\phi) \, dV \\\ I(\psi) &= - \int\_Q \psi \, c \, dS \end{aligned} \tag{10}$$

The mixed bilinear forms *C* stand for the virtual energies corresponding to the electromechanical (piezoelectric) coupling. Their definitions read:

$$\begin{aligned} \mathsf{C}(\mathfrak{v}, \mathfrak{\phi}) &= \int\_{V} \mathfrak{e}^{T}(\mathfrak{v}) \mathsf{C} \mathbf{E}(\mathfrak{\phi}) \, dV \\ \mathsf{C}(\mathfrak{\psi}, \mathfrak{u}) &= \int\_{V} \mathsf{E}^{T}(\mathfrak{\psi}) \mathsf{C}^{T} \mathfrak{e}(\mathfrak{u}) \, dV \end{aligned} \tag{11}$$

where, as before, *C* is the piezoelectric constant matrix under constant strain.

The tensorial version of the above Formulations (8)–(11) can be found in Reference [11,57]. The existence and uniqueness theorem for the above linear piezoelectricity problem can be found in Reference [60]. The proof takes advantage of the generalized Lax-Milgram lemma [61,62].

#### **3. Transition Models For Piezoelectricity**

As said in the previous sections, the transition elements under consideration serve joining the first-order shell mechanical field of a piezoelectric with the higher-order hierarchical shell models or the three-dimensional elasticity model of this field within the piezoelectric. In the electric field corresponding to the mechanical fields of both types, the same hierarchical symmetric-thickness or three-dimensional dielectric model is applied. In order to introduce the original transition models, the basic piezoelectric models that are to be joined together will be presented first.

#### *3.1. Basic Models of Linear Piezoelectricity*

The starting point for two basic models of piezoelectricity, based on the first-order and higher-order mechanical fields, is the model corresponding to the three-dimensional piezoelectricity. The matrix form of the constitutive relations of linear three-dimensional piezoelectricity, expressed by the experimentally obtainable constitutive constants, is presented in Reference [57]:

$$\begin{cases} \text{ } \mathfrak{e} = \ \mathcal{D}^{-1} \sigma - \mathfrak{e} \mathcal{E} \\\ d = -\mathfrak{c}^T \sigma + \delta \mathcal{E}. \end{cases} \tag{12}$$

Above, *<sup>D</sup>*−1, *δ*, and *c* denote the isotropic inverse matrix of elasticities (or compliance matrix), the orthotropic matrix of dielectric constants under constant stress, and the coupling orthotropic matrix of piezoelectric constants under constant stress, respectively. The six-component column vectors *σ*, *ε* represent the symmetric stress and strain tensors, *<sup>σ</sup>ij* and *<sup>ε</sup>ij*, *i*, *j* = 1, 2, 3, while *d* and *E* are the column three-component electric displacement and electric field vectors, with components *dj*, *Ej*, *j* = 1, 2, 3. The mentioned tensors and vectors are defined locally, namely in the directions consistent with the axes of piezoelectric orthotropy. In the presentation below, the third local direction is the polarization direction of the piezoelectric. In the case of the thick- or thin-walled structures, this direction coincides with the transverse one.

The standard definition of the isotropic elasticity matrix reads:


where *D* = *E*/[(1 + *ν*)(<sup>1</sup> − <sup>2</sup>*ν*)]. As it can be seen above, the non-zero terms of the matrix are expressed by Young's modulus *E* and Poisson's ratio *ν*.

The dielectricity constant matrix under constant stress reads:

$$
\mathcal{S} = \begin{bmatrix}
\delta\_{11} & 0 & 0 \\
0 & \delta\_{22} & 0 \\
0 & 0 & \delta\_{33}
\end{bmatrix} \tag{14}
$$

where *δii*, *i* = 1, 2, 3 stand for the orthotropic permittivity constants.

In the case of the piezoelectricity constant matrix under constant stress, one has

$$\mathbf{c} = \begin{bmatrix} 0 & 0 & c\_{13} \\ 0 & 0 & c\_{23} \\ 0 & 0 & c\_{33} \\ 0 & 0 & 0 \\ 0 & c\_{52} & 0 \\ c\_{61} & 0 & 0 \end{bmatrix}' \tag{15}$$

where *c*13, *c*23, *c*33, *c*52, and *c*61 are non-zero piezoelectric constants within this matrix.

3.1.1. Hierarchical or Three-Dimensional Models

For both the cases, when the transverse mechanical and electric fields are polynomially constrained through the thickness of the thin-walled piezoelectric member or the transverse direction is treated in the same way as the longitudinal directions of the three-dimensional piezoelectric member, the same constitutive relations hold.

In order to present these relations, let us convert the Equation (12) to the form suitable for the variational and finite element formulation of the general piezoelectric problem. This needs inversion of the first Equation (12) and its substitution into the second one. The resultant form reads: 

$$\begin{cases} \sigma &= \, ^\text{D}\varepsilon - \mathcal{C}E \\ \, d &= \, ^\text{C}\varepsilon + \gamma E \end{cases} \tag{16}$$

where the following definitions of the orthotropic matrix *C* of piezoelectric constants under constant strain and the isotropic matrix *γ* of dielectric constants under constant strain were introduced:

$$\begin{aligned} \gamma &= \, \, \delta - \mathfrak{c}^{\mathrm{T}} \mathrm{D} \mathfrak{c} \\ \mathfrak{C} &= \, -\mathrm{D} \mathfrak{c} .\end{aligned} \tag{17}$$

The first of the above definitions gives the following non-zero terms of the matrix *γ*:

$$
\gamma = \begin{bmatrix}
\gamma\_{11} & 0 & 0 \\
0 & \gamma\_{22} & 0 \\
0 & 0 & \gamma\_{33}
\end{bmatrix} = \begin{bmatrix}
\delta\_{11} - c\_{61}D\_{66}c\_{61} & 0 & 0 \\
0 & \delta\_{22} - c\_{52}D\_{56}c\_{52} & 0 \\
& & \delta\_{33} - c\_{13}(D\_{11}c\_{13} + D\_{12}c\_{23} + D\_{13}c\_{33}) \\
0 & 0 & -c\_{23}(D\_{21}c\_{13} + D\_{22}c\_{23} + D\_{23}c\_{33}) \\
& & -c\_{33}(D\_{31}c\_{13} + D\_{32}c\_{23} + D\_{33}c\_{33})
\end{bmatrix}.
\tag{18}
$$

The second definition leads to the following structure of the matrix *C*:

$$\mathbf{C} = \begin{bmatrix} 0 & 0 & C\_{13} \\ 0 & 0 & C\_{23} \\ 0 & 0 & C\_{33} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & C\_{52} & 0 \\ C\_{61} & 0 & 0 \end{bmatrix} = -\begin{bmatrix} 0 & 0 & D\_{11}c\_{13} + D\_{12}c\_{23} + D\_{13}c\_{33} \\ 0 & 0 & D\_{21}c\_{13} + D\_{22}c\_{23} + D\_{23}c\_{33} \\ 0 & 0 & D\_{31}c\_{13} + D\_{32}c\_{23} + D\_{33}c\_{33} \\ 0 & 0 & 0 \\ 0 & D\_{55}c\_{52} & 0 \\ D\_{66}c\_{61} & 0 & 0 \end{bmatrix} . \tag{19}$$

3.1.2. First-Order Piezoelectric Model Stress State

In the case when the mechanical field of the symmetric-thickness thick- or thin-walled member is linearly constrained through the thickness, the plane stress assumption of the first-order shell model have to be introduced into the three-dimensional piezoelectric constitutive relations. This means that the third row equation of the first constitutive relation (16), namely:

$$\sigma\_{33} = \frac{E}{(1+\nu)(1-2\nu)}[\nu\varepsilon\_{11} + \nu\varepsilon\_{22} + (1-\nu)\varepsilon\_{33}] - C\_{33}E\_{3\nu} \tag{20}$$

has to be combined with the assumption *σ*33 = 0, so as to give

$$
\varepsilon \gg = -\frac{\nu}{1-\nu} (\varepsilon\_{11} + \varepsilon\_{22}) + \frac{(1+\nu)(1-2\nu)C\_{33}}{E(1-\nu)} E\_3. \tag{21}
$$

Now, the condition (20) can be rewritten in the following alternative identity-equation form

$$\sigma\_{33} = \frac{E}{(1+\nu)(1-2\nu)} \left\{ \nu \varepsilon\_{11} + \nu \varepsilon\_{22} + (1-\nu) \left[ \frac{-\nu}{1-\nu} (\varepsilon\_{11} + \varepsilon\_{22}) + \frac{(1+\nu)(1-2\nu)\mathbb{C}\_{33}}{E(1-\nu)} E\_3 \right] \right\} - \mathbb{C}\_{33} E\_3 = 0. \tag{22}$$

Substitution of the relation (21) into the first relation (16), arithmetic manipulations on the terms containing *E*3 (with use of isotropic definitions of the non-zero terms of *Dij*, *i*, *j* = 1, 2, . . . , 6 from (13)) lead to:

$$
\begin{bmatrix}
\sigma\_{11} \\
\sigma\_{22} \\
\sigma\_{33} \\
\sigma\_{12} \\
\sigma\_{23} \\
\sigma\_{31} \\
\sigma\_{31}
\end{bmatrix} = \begin{bmatrix}
D\_{11} & D\_{12} & D\_{13} & 0 & 0 & 0 \\
D\_{21} & D\_{22} & D\_{23} & 0 & 0 & 0 \\
D\_{31} & D\_{23} & D\_{33} & 0 & 0 & 0 \\
0 & 0 & 0 & D\_{44} & 0 & 0 \\
0 & 0 & 0 & 0 & D\_{55} & 0 \\
0 & 0 & 0 & 0 & 0 & D\_{66}
\end{bmatrix} \begin{bmatrix}
\varepsilon\_{11} \\
\varepsilon\_{22} \\
\varepsilon\_{12} \\
\varepsilon\_{12} \\
\varepsilon\_{23} \\
\varepsilon\_{31}
\end{bmatrix} - \begin{bmatrix}
0 & 0 & \mathbf{C}\_{13} - \frac{\nu}{1-\nu} \mathbf{C}\_{33} \\
0 & 0 & \mathbf{C}\_{23} - \frac{\nu}{1-\nu} \mathbf{C}\_{33} \\
0 & 0 & 0 & \mathbf{0} \\
0 & 0 & 0 & \mathbf{0} \\
\mathbf{C}\_{61} & \mathbf{0} & 0
\end{bmatrix} \begin{bmatrix}
E\_{1} \\
E\_{2} \\
E\_{3}
\end{bmatrix}.
\tag{23}
$$

The alternative form of the above equation can be obtained after arithmetic manipulations on the terms containing *ε*11 and *ε*22. This leads to the change of the constitutive constants of the matrix *D*, with zero contributions of its third row and third column. This alternative form requires either retaining *ε*33 as the third term of the strain vector or removal of: the third row and the third column in *D*, the third row in *ε*, and the third row in *C*.

Taking the plane stress assumption into account in the second constitutive relation (16) needs: introduction of (21) into this relation, arithmetic manipulations on the terms containing *E*3, and summation of the contributions from the terms containing *ε*11 and *ε*22. The latter operation leads to the change of the piezoelectric constants in *<sup>C</sup>T*, including zero contributions of the third column of this matrix. All these lead to:

$$
\begin{bmatrix} d\_1 \\ d\_2 \\ d\_3 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & \mathcal{C}\_{61} \\ 0 & 0 & 0 & 0 & \mathcal{C}\_{52} & 0 \\ \mathcal{C}\_{13} - \frac{\nu}{1-\nu} \mathcal{C}\_{33} & \mathcal{C}\_{23} - \frac{\nu}{1-\nu} \mathcal{C}\_{53} & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \varepsilon\_{33} \\ \varepsilon\_{12} \\ \varepsilon\_{23} \\ \varepsilon\_{31} \end{bmatrix} + \begin{bmatrix} \gamma\_{11} & 0 & 0 \\ 0 & \gamma\_{22} & 0 \\ 0 & 0 & \gamma\_{33} + \frac{C\_{33}^2}{D(1-\nu)} \\ 0 & 0 & \gamma\_{33} + \frac{C\_{33}^2}{D(1-\nu)} \end{bmatrix} \begin{bmatrix} E\_1 \\ E\_2 \\ E\_3 \end{bmatrix}. \tag{24}
$$

Note that the third column of *C<sup>T</sup>* and the third row of *ε* can be removed from the above relation. The alternative form of (24) does not need any changes in *C<sup>T</sup>* from (16) but requires replacement of *ε*33 with −*ν* 1−*<sup>ν</sup>* (*<sup>ε</sup>*11+*ε*22), as in the Equation (23).

The chosen and written down forms (23) and (24) of the constitutive relations are suitable for further considerations in the next sections of the paper.
