3D mixed finite elements for curved, flat piezoelectric structures

The Tangential-Displacement Normal-Normal-Stress (TDNNS) method is a finite element method that was originally introduced for elastic solids and later extended to piezoelectric materials. It uses tangential components of the displacement and normal components of the normal stress vector as degrees of freedom for elasticity. For the electric field, the electric potential is used. The TDNNS method has been shown to provide elements which do not suffer from shear locking. Therefore thin structures (e.g. piezoelectric patch actuators) can be modeled efficiently. Hexahedral and prismatic elements of arbitrary polynomial order are provided in the current contribution. We show that these elements can be used to discretize curved, shell-like geometries by curved elements of high aspect ratio. The order of geometry approximation can be chosen independently from the polynomial order of the shape functions. We present two examples of curved geometries, a circular patch actor and a radially polarized piezoelectric semi-cylinder. Simulation results of the TDNNS method are compared to results gained in ABAQUS. We obtain good results for displacements and electric potential as well as for stresses, strains and electric field when using only one element in thickness direction.


Introduction
The direct and reverse piezoelectric effects allow to transform mechanical energy into electrical energy and vice versa. Moreover, modern piezoceramics are cheap in production and offer fast and accurate response. Therefore, piezoelectric materials are the method of choice in the realization of smart structures.
When it comes to design and control of piezoelectric structures, simulation results are of high interest. A popular method to get approximate solutions is the finite element (FE) method. Allik and Huges [1] and later Lerch [2] carried out first simulations using volume elements based on the principle of virtual work. These elements are an extension of standard small-strain mechanical elements by degrees of freedom for the electric potential. However, piezoelectric structures are typically flat (e.g. patch actors), and therefore in order to avoid locking effects, their discretization into well-shaped elements leads to a computational system with a huge number of unknowns -even for simple applications.
Nowadays, there are two ways to avoid this problem: the first is the derivation of equations for layered beams, plates and shells. The second is the design of locking-free volume elements. In any way, some relaxation technique is necessary to avoid shear locking. Prominent methods of relaxation are reduced integration, assumed strains, or mixed methods involving stresses, strains or dielectric displacements as additional independent unknowns. Also, a higher order of polynomial approximation is often feasible.
Mixed methods are probably the most involved of those methods mentioned above. Both from a theoretical point of view, when proving uniqueness and convergence results, as well as from an implementational point of view, when additional unknowns have to be realized in a computational code, the mixed methods require intricate treatment. However, when correctly done, mixed methods yield most accurate results, as contributions from various authors demonstrate. Among the mixed methods, we cite hybrid stress elements [3,4], a solid shell element based on a Hu-Washizu formulation [5] or the multi variable formulation [6]. Reissner-type mixed zigzag formulations were successfully used by [7,8,9]. High-order methods comprise, among many others, isogeometric and hierarchical elements from the group around Gabbert [10,11]. A high-order element with zigzag function is proposed by Polit et al. [12].
Pechstein and Schöberl [13] introduced an arbitrary order mixed FE-method for linear elasticity based on a Hellinger-Reissner formulation, which uses tangential displacements and normal normal stresses as degrees of freedom (TDNNSmethod). This method was shown to be locking-free [14], and later extended to piezoelectric materials [15]. In the current contribution, we want to show the capability of the method. We discuss high-order geometry approximation, such that curved geometries can be treated accurately by a low number of curvilin-ear elements. Moreover, we address the problem of eigenvalue computation for the mixed method. We propose to use an inverse iteration [16] for the coupled electromechanical problem. This paper is organized as follows: below, we introduce some notation on tensor products. In Section 2, the problem of linear piezolelasticity is introduced. Section 3 provides a short introduction into the variational formulation of the mixed TDNNS method that is the foundation of the current contribution. The implementation of the underlying shape functions and their derivatives is treated in Section 4. The coupled eigenvalue problem is described in Section 5, and an inverse iteration is proposed for its solution. Finally, numerical results illustrating the capability of the method are presented in Section 6.
Notation In the current contribution, we use tensor calculus. Vectorial and tensorial quantities are denoted by boldface symbols. When summing over tensor components we use Einstein's summation convention. The dot denotes contraction, while ⊗ s is the symmetric tensor product: The gradient of a scalar field a and a vector field b is denoted by the nabla operator. We indicate with respect to which coordinates the differentiation is carried out, We may omit the index, if the derivation is with respect to x. The divergence operator of a vector field is defined in the standard way, for a tensor field it is understood row-wise,

Linear piezoelasticity
We consider an elastic in part piezoelectric solid body in the domain Ω ⊂ R 3 , subjected to body forces f and free of body charges.We assume to stay in the regime of small deformations, where Voigt's theory of linear piezoelasticity [17] is appropriate. We are interested in the displacement field u and in the electric potential φ. From these quantities the electric field E = −∇φ and the linear strain tensor ε = 1 2 (∇u + ∇u T ) are derived. For the static case the mechanical and electrical balance equations in differential form are given by with the Cauchy stress tensor σ and the dielectric displacement D. The mechanical boundary conditions are and the electrical boundary conditions We are concerned with the most common variant of constitutive equations, which are often referred to as Voigt's linear theory of piezoelectricity [17]. In our notation, they read with C E the elasticity tensor measured under constant electric field, ǫ ε the dielectric tensor at constant strain and the piezoelectric permittivity tensor e. We use the following alternative variant of the material law Of course the material parameters are connected, their relations are

TDNNS for piezoelasticity
In this section we briefly introduce the "Tangential Displacement Normal Normal Stress" (TDNNS) finite element method. The TDNNS finite element method was introduced by Pechstein and Schöberl [13,18] for elasticity and later extended to piezoelasticity [15]. It is a mixed finite element method based on Reissner's principle which considers displacements and stresses as unknowns. A more detailed description can be found in [15]. We assume T = {T } to be a finite element mesh of the domain Ω. The unit vector n denotes the outer normal on element or domain boundaries ∂T or ∂Ω. In general on a (boundary) surface a vector field v can be split into a normal component v n = v · n and a tangential component v t = v − v n n. A tensor field τ has a normal vector on a surface τ n = τ · n, which analogously can be split into a normal component τ nn and tangential component τ nt .

Finite element formulation for elastic problem
The TDNNS finite element method uses the tangential component of the displacement u t and the normal component of the stress vector σ nn as degrees of freedom. Those quantities are continuous across the element interfaces. Note that the normal displacement u n is discontinuous. For the displacements, tangentialcontinuous elements introduced by Nédélec [19] are used. Elements for the stresses are normal-normal-continuous and were introduced by Pechstein and Schöberl [13]. For the electric potential, standard continuous (e.g. nodal or hierarchical) finite elements are used. We assume the tangential continuous Nedelec elements for the displacements and continuous nodal or hierarchical elements for the electric potential on hybrid meshes to be well known. The definition of arbitrary-order prismatic or hexahedral normal-normal continuous stress elements is postponed to Section 4. Other element types were introduced in [13]. At present, it is only necessary to assume that admissible stresses σ are such that σ nn is continuous across all element interfaces.
The variational formulation of the elastic problem based on Reissner's principle is the following: find u satisfying u t = 0 on Γ 1 and σ satisfying σ nn = 0 on Γ 2 such that for all admissible virtual displacements δu and virtual stresses δσ which satisfy the corresponding homogeneous essential boundary conditions. In (13) instead of the integral Ω ε(u) : σ dΩ, the duality product ε(u), σ occurs. This is due to the fact that the displacement field u is discontinuous, and therefore gaps between elements may arise. On finite element meshes, the duality product can be computed element wise by surface and volume integrals . Due to (additional) distributional parts of the strain, the surface integrals in (14)- (15) do not vanish. Only if the stress field is normal-normal continuous, the duality product is well defined Equations (14)- (15) can be shown by partial integration, using continuity of normal normal stress and tangential displacement and is shown in detail in [13] and [18].

Finite element formulation for piezoelastic problem
To get equations for the (fully) coupled problem of piezoelasticity we eliminate the dielectric displacement in the balance equations (4)-(5) and use the constitutive equations in (10)-(11) to get We multiply (16) by a virtual stress δσ, (17) by a virtual displacement δu and (18) by a virtual electric potential δφ to get a variational formulation. The virtual quantities have to satisfy the corresponding homogeneous boundary conditions δu t = 0 on Γ 1 , δσ nn = 0 on Γ 2 and δφ = 0 on Γ 3 . After integration using the identities in (14) -(15) we get from (16) and (17) Integration by parts in (18), taking into account the homogeneous boundary conditions Γ 3 (δφ = 0) and the natural one on Γ 4 (no surface charges) leads to Summing up we finally arrive at The performance of thin prismatic elements has been shown in [15]. They have been shown to be free from locking when using at least a order of k = 2 for the electric potential and linear elements for the mechanical quantities (stress and displacement).

Curvilinear elements
To implement the equations from the previous section in a computational code, we use the well known concept of reference elements. The reference element is of unit size and not distorted. Contrarily, elements in the finite element mesh are in general distorted, and may be curved in order to enhance geometry approximation. Integration of virtual work is then carried out after transformation to the reference element. Also shape functions are defined for the reference element.
In Section 4.1, we define the reference element for elements of prismatic and hexahedral shapes. Then we introduce shape functions for the stresses on these elements. Note that, for prismatic elements, we use less shape functions than proposed in the original paper [14], while shape functions for hexahedral elements are introduced for the first time. In Section 4.2, we describe in detail how quantities on the reference element are transformed to a distorted element in the mesh. Displacement shape functions are transformed by a covariant transformation, which preserves tangential degrees of freedom. From this transformation, we derive the correct transformation for the strain on a curvilinear element. Moreover, we present the correct transformation for the stress shape functions, as well as the transformation of the divergence of the stresses.

Reference element and shape functions
In the following, we will provide shape functions of arbitrary order for displacements, stresses and electric potential on the reference hexahedral and prism. As we use the tensor product structure of these elements, we first describe the reference segment and reference triangle. All quantities associated to reference elements are denoted by a hat. The shape functions are then transformed by according transformations described in Section 4.2 to some possibly distorted element in the mesh.
For the electric potential, we use standard continuous hierarchical elements. For the displacement elements, we use tangential continuous elements well-known from electromagentics [19]. All elements are implemented in the open-source software package Netgen/NGSolve [20]. The exact basis function implemented there are described in detail in [21]. We keep close to the notation adopted in this latter reference. The stress basis functions were introduced for tetrahedra in [13] and for prismatic elements in [14]. In the following, we provide an adapted set of stress basis functions for prisms and stress basis functions for hexahedral elements.
The reference segment We use the unit reference segmentT seg = [0, 1]. For reference coordinatex, we define the barycentric coordinates or "hat basis func- Figure 1: The various unit elements.
and the family of Legendre polynomials,q i =l i (x), where i indicates the polynomial order.
The reference triangle The unit triangle shall be denoted byT trig = {(x,ŷ) : x ≥ 0,ŷ ≥ 0,x +ŷ ≤ 1}. For an enumeration of vertices and edges we refer to Figure 1. Also on the triangle, we define barycentric coordinateŝ We need a family of polynomials on the edge E γ between vertices V α and V β , that is extended to the whole triangle, e.g. the family of scaled Legendre polynomialŝ A family of bivariate polynomials on the triangle of polynomial order i + j, can e.g. be realized byq Of course any other family of polynomials on the triangle may be used. The conditioning of the finite element matrices for high polynomial orders is affected by this choice, but this is not within the scope of the current contribution.
On the triangle, we introduce three unit tensor fields, that have unit normal stress on one edge and zero normal stress on the other edges. Multiplying these unit tensors by (scalar) polynomials, we arrived at the desired hierarchical shape functions on the triangle, they can be found in [13]. We set the unit tensorŝ The reference prism The reference prism is constructed as a tensor product of the reference triangle (forx andŷ coordinates) and the reference segment (for z coordinate): The reference prism is depicted in Figure 1.
We define the local reference shape functions for the element of order p. On the prism, we have shape functions associated to each of the three quadrilateral faceŝ Fxŷ 1 ,Fxŷ 2 andFxŷ 3 , as well as shape functions associated to the two triangular faceŝ Fẑ 1 andFẑ 2 , and shape functions that have zero normal stress and are associated with the element itself. The last class of shape functions is element-local and can be eliminated directly at assembly. We use • for each triangular faceFẑ γ , 1 2 (p + 1)(p + 2) shape functionŝ • three types of element-local shape functions with vanishing normal stress, each associated to one block of the symmetric three-by-three tensor, Note that some components of the shape functions are of order p + 1, however the element normal stresses are of order p.
The reference hexahedron The reference hexahedron is again constructed in tensor product style:T hex =T seg ⊗T seg ⊗T seg . The reference hexahedron is depicted in Figure 1. Again, we distinguish between shape functions that are associated with one of the quadrilateral faces, with zero normal stress on all other faces, and shape functions that have zero normal stress on all faces. Additionally, the shape functions are built in such a way that the normal stress of the face-based shape functions corresponds to the normal stress distribution of face-based shape functions for prisms. This way, it is possible to use prismatic and hexahedral elements in the same mesh. Of course, also the shape functions for tetrahedral elements are implemented in NGSolve, which are not discussed here, match accordingly, such that hybrid meshes can be used. We use the shape functions • for each of the facesFξ γ withξ ∈ {x,ŷ,ẑ} and γ = 1, 2 • the element-local shape functions with vanishing normal stresŝ

Transformations
Let nowT denote the reference element of any type (triangular, quadrilateral, prismatic, tetrahedral or hexahedral), and let T be a corresponding element in the finite element mesh. We assume Φ T (x) to be a smooth one to one mapping fromT to T . A point in the reference elementx ∈T is mapped by Φ T to some point x ∈ T . In Figure 2 this transformation is illustrated. Note, that in general Φ T is nonlinear. The Jacobian of the transformation (similar to the deformation gradient tensor) is denoted by The determinant of the Jacobian is denoted J T = det(F T ). Furthermore the Hessian of the i th component of Φ T , H i T , is given by As already mentioned, the degrees of freedom of the method are u t and σ nn . These degrees of freedom need to be preserved when shape functions are transferred to the curvilinear element. AssumeN σ andN u to be one specific displacement and stress shape function, respectively. To calculate a finite element function, we transform these shape functions like The displacement and stress field are then weighted sums of these basis functions. The first transformation for the displacement shape functions (39) is well known from finite elements for the electric field in Maxwell's equations, see e.g. [22]. The linear strain ε of a finite element function is the the weighted sum of strains of such basis functions ε(N u ). On curvilinear elements, this strain of a basis function has to be evaluated using the chain rule. In contrast to standard elements, it depends not only on the reference strainε(N u ) = 1 2 ∇xN u + ∇xN u,T but also on the shape functionN u The second part of (41) only vanishes if F T is constant, which is the case only for affine linear elements.
The transformation for stress shape functions (40) can be found in [13]. Note that compared to the Piola transformation, we have a factor of 1 For the (mechanical) balance equation (4) the divergence of the stress tensor has to be derived. Again, it consists of a weighted sum of div x (N σ ). From basic calculus for the standard Piola transformation, we deduce whereF T := 1 J T F T . By application of the product rule we get the i th component of the divergence of the stress tensor as Assuming that the shape functionsN σ as well as their divergence divxN σ are known analytically on the reference element, formula (43) allows to evaluate div x N σ for a curvilinear mesh element.

Eigenvalue problem
A wide range of technical applications involving piezoelectric structures is operated in the regime of (mechanical) vibrations e.g. noise and vibration control and damping, shape control or energy harvesting. To design and study these applications the computation of their eigenfrequencies and -forms may be required. The finite element discretization leads to a system with a huge number of eigenvalues of which typically only a few of the smallest are of interest. In our contribution we propose to use the inverse iteration procedure which allows the effective computation of a certain number of eigenvalues.

Inverse iteration
The inverse iteration procedure [23, 16,24] is designed to find the smallest eigenvalues and -vectors of the generalized eigenvalue problem with symmetric positive definite stiffness matrix A and mass matrix M. For a given λ n k and q n k , the next iterate q n+1 k is computed by solving a linear system An approximation of the eigenvalue is given by the Rayleigh quotient For a suitable start vector q 0 , which is not perpendicular to the smallest eigenvector, the inverse iteration procedure converges to the smallest eigenvalue λ 1 and the corresponding Eigenvector q 1 .

Structure of assembled system
The TDNNS (mixed) method for piezoelasticity with the variational formulation (22-23) leads to a certain structure of the global assembled system. The global stiffness matrix A global is indefinite, the global mass matrix M global is positive semidefinite. The eigenvalue problem is of the form The global vector of degrees of freedom is denoted by q, where q u , q σ and q ϕ denote the degrees of freedom of displacement, stress and electric potential, respectively. The symmetric positive definite matrices C and E correspond to the mechanical compliance and electrical permittivity of the constitutive equations.
In the variational formulation these parts are represented by Ω σ : S E : δσdΩ and Ω ∇φ·ǫ·δ∇φdΩ. The matrix D corresponds the electromechanical coupling terms Ω (d : σ · ∇δφ) dΩ, while B represents ε(u), δσ , defined in (14)- (15). Note that in (47) the right hand side is 0 for the stresses and the electric potential. Therefore the global mass matrix is positive semidefinite.
In the following, we motivate why the inverse iteration procedure can be used although neither the mass nor the stiffness matrix are positive definite. We emphasize that these considerations are purely theoretical, while the procedure is then directly applied to the indefinite or semidefinite matrices in implementations. From the last line of (47) the potential degrees of freedom can be computed as Inverting (48) into the second line of (47) we get for the stresses Finally applying (49) to the fist line of (47) we get the eigenvalue problem For typical technical piezoelectric materials the matrixC can be assumed to be positive definite as it represents the compliance at constant dielectric displacement. Therefore, (49) is formally equivalent to a generalized eigenvalue problem with symmetric positive definite matrices, and can be treated by the inverse iteration procedure.
Note that the elimination of the degrees of freedom of stress an potential in (48)-(50) is only of theoretical interest, to show that the eigenvalue problem can be solved by the inverse iteration procedure. In the practical implementation the global stiffness and mass matrix A global and M global are used. An exemplary implementation of the inverse iteration procedure Netgen/NGSolve can be found in the tutorial section in [20].

Numerical Results
In the sequel we study two numerical examples. First we show the discretization of a circular piezoelectric patch actor applied to a square aluminium plate. We evaluate the displacement of the plate and the stress field for an applied static voltage as well as the eigenfrequencies and eigenforms of the assembly. Our second example is a radially polarized piezoelectric half cylinder, studied in [25]. We calculate (static) displacements and present a convergence study. For both examples we compare our results to results gained in the commercial software tool ABAQUS 6.14 [26].

Circular patch actor
Our first numerical example is a circular piezoelectric patch applied on a square plate. A schematic sketch is shown in Figure 3. The piezoelectric patch has a diameter of 15 mm and a thickness of 0.5 mm. The patch material is considered to be PTZ-H5, polarized in thickness direction. The material properties are taken from [25] and summarized in Table 1. The square aluminium plate (Young's modulus 65 GPa, Poisson ratio 0.3, density 2.7 × 10 −9 kg/mm 3 ) has a length of 25 mm and a thickness of 1 mm. One side of the plate is clamped. A potential difference ∆φ = 100 V is applied to the electrodes of the piezoelectric patch, which are located at the top and bottom of the patch actor.
PSfrag replacements For our reference calculations in ABAQUS we use a mesh with an element size 0.25 mm, consisting of four elements across the thickness of the aluminium and two elements across the thickness of the patch. We use ABAQUS 6.14 C3D20R elements (20 nodes quadratic brick elements, reduced integration) for the linear elastic aluminium plate and C3D20RE elements (20 node quadratic piezoelectric brick, reduced integration) for the piezoelectric patch. Our ABAQUS model requires approximately 830k degrees of freedom.
The mesh used for the calculation in Netgen/NGSolve is shown in Figure 3b. It consists of a very coarse prismatic mesh using only one element in thickness direction, for both the plate and the patch. On the circumference of the patch the mesh is refined by two ring layers of hexahedral elements, which have an aspect ration of ≈ 25 : 2 : 1. We use shape functions of order p = 2 for displacements and stresses and order p φ = 3 for the electric potential. Summing up, we use 12439 degrees of freedom for calculations in Netgen/NGSolve.
Note that in the TDNNS-method, contrary to elements using nodal degrees of freedom, a reduction of the thickness of the patch (or the plate) has no effect on the number of used degrees of freedom. Actually we have chosen a rather thick patchactor, in order to get accurate reference solutions at reasonable computational costs.
We compare resulting stress fields. In Figure 4 show contour plots of the stress component S 11 . In Figure 4 the stress fields calculated in Netgen/NGSolve and reference solution are shown. We as well evaluate the resulting vertical displacements of the corner point marked in Figure 3(b). While the reference displacement is −6.07325 µm, our result is −6.03698 µm. The relative difference of the results is −0.597%.  Table 2. For the short cut configuration reference eigenfrequencies were calculated in ABAQUS. The relative frequency error ∆f SC = f SC,ref Table 2. The lowest three eigenforms are plotted in Figure 5.

Radially polarized semicylinder
Our second example is a radially polarized piezoelectric semicylinder that was first presented in [25]. The semicylinder is of diameter 15 mm, thickness 1 mm and length (in z direction) 5 mm. A sketch of the geometry is shown in Figure 6. We re-use the material parameters from Table 1 in the local frame of the polarization direction. To the electrodes, located at the inner and outer cylinder surfaces, we apply a potential difference of 100 V . We compute a reference solution in ABAQUS 6.14 using 10540 C3D20Eelements (four across the thickness, 17 across the length and 155 across the circumference), in total 210k degrees of freedom. The reference tip displacement at the marked corner point in Figure 6(b) is u tip ref = 2.75525 µm. The mesh used for calculations in Netgen/NGSolve is shown in Figure 6b. Again a very coarse prismatic mesh consisting of one element along the length of semicylinder and 10 elements along the circumference is used. This corresponds to a mesh size h = 5. A second (finer) mesh with two elements in z-direction For both meshes we use only one element in thickness direction. We performed calculations with different orders of interpolation k for displacements and stresses, interpolation order for the electric potential is k + 1. For each calculation the tip displacement u tip at the corner point is evaluated. We evaluate the relative difference of results defined as δ rel = u tip /u tip ref − 1. The results of this convergence study are shown in Table 3 together with the number of degrees of freedom n dof used in Netgen/NGSolve. We see that the accuracy of the solution increases both, when decreasing the mesh size or increasing the polynomial order. More detailed convergence studies can be found in [15].

Conclusion
In our contribution we have provided prismatic and hexahedral elements of arbitrary polynomial order for the discretization of flat possibly curved piezoelectric geometries by the TDNNS-method. The inverse iteration has been proposed to calculate eigenvalues and -forms. In our numerical examples we have shown, that accurate simulation of flat curved piezoelectric structures can be done with only one element in thickness direction. Appropriate results were achieved for displacements and stresses as well as for frequencies, while the number of degrees of freedom could be reduced significantly.