Exact geometry SaS solid-shell element for 3D stress analysis of FGM piezoelectric structures

Abstract A hybrid-mixed functionally graded material (FGM) piezoelectric four-node solid-shell element through the sampling surfaces (SaS) method is proposed. The SaS formulation is based on choosing inside the shell N SaS parallel to the middle surface in order to introduce the displacements and electric potentials of these surfaces as fundamental shell unknowns. Such choice of unknowns with the use of Lagrange polynomials of degree N-1 in through-thickness interpolations of the displacements, strains, electric potential, electric field and material properties leads to a robust FGM piezoelectric shell formulation. The inner SaS are located at Chebyshev polynomial nodes that make it possible to minimize uniformly the error due to Lagrange interpolation. To implement the effective analytical integration throughout the element, the extended assumed natural strain (ANS) method is employed. As a result, the piezoelectric four-node solid-shell element exhibits a superior performance in the case of coarse meshes. To circumvent shear and membrane locking, the hybrid stress-strain solid-shell formulation via the Hu-Washizu variational principle is employed. The developed solid-shell element could be useful for the 3D stress analysis of FGMstructures because the SaS method allows obtaining the solutions with a prescribed accuracy, which asymptotically approach the exact solutions of electroelasticity as the number of SaS tends to infinity.


Introduction
In recent years, a considerable work has been carried out on continuum-based finite elements [1][2][3][4][5] that can handle the analysis of shells satisfactorily. These elements are defined by two layers of nodes at the bottom and top surfaces with three translational degrees of freedom (DOF) per node and known as 6-parameter solid-shell elements. They were applied effectively for the analysis of coupled electro-mechanical systems [6][7][8][9][10][11][12][13][14]. However, the 6-parameter solid-shell element formulation based on the complete constitutive equations of piezoelectricity is deficient because thickness locking occurs. This is due to the fact that the linear displacement field in the thickness direction results in a constant transverse normal strain, which causes artificial stiffening of the shell element in the case of non-vanishing Poisson's ratios. To prevent thickness locking, the 3D constitutive equations can be modified using the generalized plane stress conditions [7,12]. The hybrid stress method [6,14] in which the transverse normal stress is assumed to be constant through the shell thickness and the most popular enhanced assumed strain (EAS) method in which the transverse normal strain is enriched in the thickness direction by a linear term [2] are also utilized [8][9][10][11]13].
An effective way of using the 3D constitutive equations is to employ the solid-shell model with seven translational DOF [15][16][17][18][19][20][21]. The 7-parameter shell formulation is based on choosing six displacements of the bottom and top surfaces and the transverse displacement of the middle surface as basic shell unknowns. Such formulation is optimal with respect to the number of DOF. To circumvent locking phenomena, the ANS method [15][16][17][18][19], EAS method [17] and hybrid-mixed method [20,21] were employed.
The more general 9-parameter solid-shell formulation utilizes nine displacements of external and middle surfaces as shell unknowns [22]. Such choice of unknowns with the consequent use of Lagrange polynomials of the second order in through-thickness approximations of the displacements and strains leads to a robust higher-order shell formulation. Moreover, this model allows the derivation of objective strain-displacement equations, which ex-actly represent rigid-body motions of the shell in any convected curvilinear coordinate system. Taking into account that the displacement vectors of reference surfaces are resolved in the middle surface frame, the higher-order shell formulation with nine translational DOF is very promising for developing the exact geometry or geometrically exact (GeX) solid-shell elements. The term GeX means that the parametrization of the middle surface is known a priori and, therefore, the coefficients of the first and second fundamental forms and Christoffel symbols are taken exactly at element nodes. The 9-parameter solid-shell element has been developed in contributions [23,24] to solve the coupled problem of electroelasticity for piezoelectric structures.
It is well known that aforementioned isoparametric and GeX piezoelectric solid-shell elements [6-14, 20, 21, 23, 24] do not describe properly transverse components of the stress tensor and electric displacement. To calculate them, the higher-order shell models have to be employed. The robust GeX higher-order piezoelectric shell elements with the various number of DOF per node are proposed by Carrera and his coauthors [25][26][27] through Carrera's unified formulation. These finite elements are applied to the 3D electro-mechanical analysis of piezoelectric plates and cylindrical shells. However, the authors [26] report that the piezoelectric nine-node element based on the fourth-order shell theory accounting for thickness stretching does not provide fulfilling the boundary conditions for the transverse normal stress and electric displacement on the bottom and top surfaces especially in the case of thin shells.
Recently, this problem has been solved by means of the GeX four-node piezoelectric solid-shell element [28] using a SaS method. The SaS piezoelectric shell formulation [29,30] is based on choosing inside the shell N not equally spaced SaS Ω 1 , Ω 2 ,. . . , Ω N parallel to the middle surface and located at Chebyshev polynomial nodes [31] in order to introduce the displacements and electric potentials of these surfaces as fundamental shell unknowns, where N ≥ 3. Such choice of unknowns with the use of Lagrange polynomials of degree N − 1 in assumed distributions of the displacements, strains, electric potential and electric field through the thickness yields a very compact piezoelectric shell formulation.
The functionally graded materials (FGMs) are a new class of advanced materials in which the material properties vary continuously from point to point. This property is achieved by varying the volume fraction of constituents. One of the advantages of the monotonic variation of the volume fraction of constituent phases is the elimination of stress discontinuity that is often encountered in laminated composites. Nowadays, the solid-shell elements are widely used in bending, buckling and vibration analyses of FGM shell structures because of their advantages compared with other numerical techniques. The first contribution to this area is an isoparametric FGM four-node solid-shell element proposed by El-Abbasi and Meguid [32] that generalizes the ANS solid-shell element with seven translational DOF per node [18]. Arciniega and Reddy [33] developed the GeX FGM solid-shell element based on the 7-parameter shell formulation accounting for thickness stretching. By increasing the p-level of the finite element approximation within each element, they eliminate shear and membrane locking. This technique was utilized by Payette and Reddy [34] in a 7-parameter spectral/hp solid-shell element formulation. The efficient GeX nine-node FGM shell elements based on the second-, third-and fourth-order models via Carrera's unified formulation have been developed in [35,36]. To circumvent shear and membrane locking, the ANS concept is employed. The most general finite element on the basis of the fourthorder theory [35] exhibits a high accuracy for transverse stresses except for the transverse normal stress especially in thin shell limits. To solve a problem, the FGM quadrilateral solid-shell element through the SaS concept has been developed [37]. This finite element extends authors' SaS solid-shell elements [38,39] to FGMs.
The purpose of this paper is to combine the FGM quadrilateral element [37] with the GeX piezoelectric solidshell element [28]. Thus, we deal with a general GeX four-node solid-shell element with the inner SaS located at Chebyshev polynomial nodes for the coupled electromechanical analysis of FGM piezoelectric doubly-curved shells. To overcome shear and membrane locking, the assumed interpolations of stresses and displacementindependent strains are utilized through the Hu-Washizu variational principle. Such approach exhibits an excellent performance in the case of coarse mesh configurations and has computational advantages compared to conventional isoparametric hybrid-mixed solid-shell element formulations because it reduces the computational cost of numerical integration in the evaluation of the element stiffness matrix. This is due to the fact that all element matrices require only direct substitutions, i.e., no expensive numerical matrix inversion is needed. It is impossible in the framework of the isoparametric hybrid-mixed piezoelectric shell element formulation [6,7,10,11,14]. Second, the GeX piezoelectric four-node solid-shell element formulation is based on efficient analytical integration throughout the element using the extended ANS method [22,24]. The latter has a great meaning for the numerical modelling of doubly-curved shells.

SAS formulation for piezoelectric shell
Consider a piezoelectric shell of the thickness h. Let the middle surface Ω be described by orthogonal curvilinear coordinates θ 1 and θ 2 , which are referred to the lines of principal curvatures of its surface. The coordinate θ 3 is oriented along the unit vector e 3 (θ 1 , θ 2 ) normal to the middle surface (see Figure 1). Introduce the following notations: eα(θ 1 , θ 2 ) are the orthonormal base vectors of the middle surface; Aα(θ 1 , θ 2 ) are the coefficients of the first fundamental form; kα(θ 1 , θ 2 ) are the principal curvatures of the middle surface; cα = 1 + kα θ 3 are the components of the shifter tensor; c I α (θ 1 , θ 2 ) are the components of the shifter tensor at SaS Ω I given by where θ I 3 are the transverse coordinates of SaS defined as where N is the number of SaS. Here and in the following developments, the indices I, J, K identify the belonging of any quantity to the SaS and run from 1 to N; Latin indices i, j, k, l range from 1 to 3; Greek indices α, β range from 1 to 2. It is worth noting that the inner SaS are located at the roots of the Chebyshev polynomial of degree N − 2 [40]. The through-thickness SaS approximations [29] can be written in a compact form as where u i , ε ij , σ ij , ϕ, E i are the displacements, strains, stresses, electric potential and electric field; u I i (θ 1 , θ 2 ), ε I ij (θ 1 , θ 2 ), σ I ij (θ 1 , θ 2 ), ϕ I (θ 1 , θ 2 ), E I i (θ 1 , θ 2 ) are the displacements, strains, stresses, electric potential and electric field of SaS Ω I ; L I (θ 3 ) are the Lagrange basis polynomials of degree N − 1 defined as , In the orthonormal basis e i , the relations between strains and displacements of SaS [31] are expressed as where λ I iα (θ1, θ 2) are the strain parameters of SaS; β I i (θ1, θ 2) are the values of the derivative of displacements with respect to thickness coordinate on SaS: where the symbol (. . .) ,i stands for the partial derivatives with respect to coordinates θ i ; M J = L J ,3 are the polynomials of degree N − 2; their values on SaS are calculated as It is seen from (8) that the key functions β I i of the SaS piezoelectric shell formulation are represented as a linear combination of displacements of SaS u J i . In the orthonormal basis e i , the relations between the electric field and electric potentials of SaS [29] are written as One can see that the normal component of the electric field E I 3 is represented as a linear combination of electric potentials of SaS ϕ J that is similar to (8).

Hu-Washizu variational principle for FGM piezoelectric shell
The proposed hybrid-mixed FGM piezoelectric solid-shell element is based on the modified Hu-Washizu variational principle of piezoelectricity in which displacements, strains, stresses and electric potential are utilized as independent variables [28]: is the infinitesimal volume element; ε ij and η ij are the displacement-dependent and displacement-independent strains; C ijkl , e kij and ∈ ij are the elastic, piezoelectric and dielectric constants; W is the work done by external electro-mechanical loads. As usual, the summation on repeated Latin indices is implied. Following the SaS technique (3), we introduce the next assumption of the hybrid-mixed solid-shell element formulation. Assume that the displacement-independent strains are distributed through the thickness of the shell by where η I ij (θ 1 , θ 2 ) are the displacement-independent strains of SaS.
Finally, we accept the last assumptions of the SaS FGM piezoelectric shell formulation. Let the material constants be distributed through the shell thickness as follows: where C I ijkl (θ 1 , θ 2 ), e I kij (θ 1 , θ 2 ) and ∈ I ij (θ 1 , θ 2 ) are the values of elastic, piezoelectric and dielectric constants on SaS.
Substituting the through-thickness distributions (3), (13)- (16) in (12) and introducing the weighted coefficients one can write the Hu-Washizu functional in terms of only SaS variables as

Hybrid stress-strain FGM piezoelectric solid-shell element formulation
The finite element formulation is based on a simple interpolation of the shell via GeX four-node piezoelectric solidshell elements where Nr (ξ1, ξ 2) are the bilinear shape functions of the element; u I ir and ϕ I r are the displacements and electric potentials of SaS Ω I at element nodes; ξ 1 , ξ 2 are the normalized curvilinear coordinates θ 1 , θ 2 ( Figure 2); the nodal index r runs from 1 to 4.
To implement the efficient analytical integration throughout the finite element, the extended ANS method [22,24] can be adopted to interpolate the displacement-dependent strains and the electric field where ε I ijr and E I ir are the strains and electric field of SaS at element nodes. Remark 1. The main idea of such approach can be traced back to the ANS method developed by many scientists [41][42][43][44] to cure the isoparametric finite elements from shear and membrane locking. In contrast to the conventional formulation, we treat the term ANS in a broader sense. In the GeX four-node solid-shell element formulation, the displacement-dependent strains of SaS are assumed to vary bilinearly throughout the biunit square in (ξ 1 , ξ 2 )space. The extended ANS method (22) and (23) makes it possible to utilize the element nodes as sampling points that helps to avoid the use of Gauss numerical integration.

Remark 2.
In order to circumvent curvature thickness locking for the isoparametric four-node solid-shell element, Betsch and Stein [45] proposed to utilize the bilinear interpolation (22) for the transverse normal strain. It is apparent that curvature thickness locking is not related to the GeX four-node solid-shell element because it can handle the arbitrary geometry of surfaces properly. We advocate the use of the extended ANS method for all components of the strain tensor to implement the effective analytical integration throughout the element.
The strain vectors of SaS at element nodes can be expressed as where B I u r are the constant inside the element matrices of order 6 × 12N; U is the element displacement vector given by The electric field vectors of SaS at element nodes are where B I ϕ r are the constant inside the element matrices of order 3 × 4N; Φ is the element electric potential vector defined as For further developments, it is convenient to write the ANS interpolation (22) in the following form: where Here, and below the indices r 1 and r 2 run from 0 to 1.
The same improvement has to be done concerning the ANS interpolation (23), that is, where To overcome shear and membrane locking and obtain no spurious zero energy modes, the robust 12-parameter stress interpolation [28] is utilized The similar interpolation is adopted for displacementindependent strains: Substituting interpolations (20), (21), (28), (30), (32) and (33) in the Hu-Washizu variational principle (11) and (18), replacing the metric product A 1 A 2 in surface integrals by its value at the element center and integrating analytically throughout the finite element, the following equilibrium equations of the GeX hybrid-mixed FGM piezoelectric four-node solid-shell element are obtained: for r 1 + r 2 < 2, where Fu and F ϕ are the mechanical and electric surface vectors. Eliminating vectors η I r1r2 and σ I r1r2 from Eqs. (35)-(38) and accounting for that det(Γ IJ ) ≠ 0, we arrive at the system of linear equations where Kuu, K uϕ , K ϕu = K T uϕ and K ϕϕ are the mechanical, piezoelectric and dielectric stiffness matrices defined as

Remark 3.
As pointed out earlier, all element matrices are evaluated with no expensive numerical matrix inversion that is impossible in conventional isoparametric hybridmixed piezoelectric solid-shell element formulations.

Numerical examples
The performance of the proposed GeX four-node FGM piezoelectric solid-shell element on the basis of the SaS technique denoted by the GeXSaS4 element is assessed through the use of exact and numerical solutions of piezoelectricity extracted from the literature.

FGM piezoelectric rectangular plate under sinusoidal loading
Consider a simply supported FGM piezoelectric rectangular plate subjected to mechanical loading acting on its top surface with boundary conditions where a and b are the plate dimensions. Let the FGM constants be distributed in the thickness direction according to the exponential law [46]: where C − ijkl , e − kij and ∈ − ij are the values of elastic, piezoelectric and dielectric constants on the bottom surface, which are considered to be the same as those of the PZT-4 given in [46] and Table 1; α is the material gradient index, which can be determined as where C + ijkl , e + kij and ∈ + ij are the values of elastic, piezoelectric and dielectric constants on the top surface.
To compare the derived results with the exact solution [46] of piezoelectricity, we accept a = b = 1 m and p 0 = 1 Pa, and introduce the scaled basic variables as follows:   Due to symmetry of the problem, only one quarter of the plate is modeled by uniform meshes depicted in Figure 3. The data listed in Tables 2 and 3 show that the GeXSaS4 element makes it possible to reproduce the exact solution of piezoelectricity [46] for the FGM piezoceramic  Table 3: Results for a FGM piezoelectric square plate with a/h = 10 and α = −1 using a uniform 64 × 64 mesh.
Formulation  [47] by choosing seven SaS inside the plate is given. As can be seen, the GeXSaS4 element provides already four right digits for most basic variables since seven SaS. Figure 4 displays distributions of the transverse displacement, stresses, electric potential and electric displacement (44) through the thickness of the plate for different slenderness ratios using seven SaS and the same fine mesh. These results demonstrate the high potential of the GeXSaS4 element because the boundary conditions on bottom and top surfaces for the transverse stresses and electric displacement are satisfied correctly. Figure 5 shows the results of the convergence study due to mesh refinement through the transverse displacement, electric potential and stresses for different slenderness ratios choosing seven SaS inside the plate body. The reference values are provided by authors' exact SaS solution [47] using the same number of SaS. It is seen that the GeXSaS4 element behaves well in the case of coarse mesh configurations.

FGM piezoelectric spherical shell under electric loading
Next, we study a FGM piezoelectric spherical shell subjected to a constant potential of ϕ 0 = 1 V applied on the outer surface and consider the following boundary conditions: This problem is a good benchmark to test the proposed analytical integration schemes because in the literature there is Heyliger-Wu's exact solution [48] for the homogeneous piezoceramic spherical shell. It is supposed that the FGM constants are distributed through the thickness according to the exponential law (42), where C − ijkl , e − kij and ∈ − ij are the values of elastic, piezoelectric and dielectric constants on the inner surface, which are taken to be the same as those of the PZT-5 given in [48] and Table 1.
Owing to symmetry, one sixteenth of the shell is discretized by regular 4n × 1 meshes shown in Figure 6. The coefficients of the first and second fundamental forms and Christoffel symbols of this part of the spherical surface with a hole at the top are defined as where R = 0.045 m and θ 0 = π/18000. It is convenient to introduce the following scaled variables as functions of the dimensionless thickness coordinate: σ 11 (z) = 10 −2 × σ 11 (P, z), σ 33 (z) = 10 −2 × σ 33 (P, z), ϕ(z) = ϕ(P, z), D 3 (z) = 10 7 × D 3 (P, z), z = θ 3 /h, where P(π/4, 0) is the point belonging to a middle surface. Tables 4 and 5 list the results of the convergence study due to increasing the number of SaS inside the shell body by using a regular 64 × 1 mesh. A comparison with the     Formulation   Heyliger-Wu's solution [48] is also given. As it turned out, the GeXSaS4 element provides already four right digits for all basic variables starting respectively from nine and five SaS for thick and moderately thick spherical shells. Figures 7-9 show the through-thickness distributions of the displacement, stress, electric potential and electric displacement for different values of the slenderness ratio R/h and material gradient index α choosing nine SaS and 64×1 mesh. One can see that boundary conditions for the transverse normal stress are satisfied again properly with a high accuracy. The results of the convergence study due to mesh refinement through the normalized displacement, stress, electric potential and electric displacement for the same number of SaS are presented in Figure 10. The reference values have been obtained taking nine SaS and a fine 128 × 1 mesh. It is seen that the GeXSaS4 element behaves practically insensitive with respect to coarse meshes in thick and thin shell limits.

FGM piezoelectric cylindrical shell under opposite line loads
Finally, we study a FGM piezoelectric cylindrical shell subjected to the uniformly distributed line load p + 3 = 656.17 N/m acting on the top surface Ω + along two diametrically opposite lines as depicted in Figure 11. The piezoelectric is polarized in the thickness direction and the electrodes on bottom and top surfaces are grounded.
It is assumed that the FGM constants are distributed through the shell thickness according to the exponential law (42), where C − ijkl , e − kij and ∈ − ij are the material constants on the bottom surface, which are considered to be the same as those of the PZT-4 given in [49,50] and Table 1. The geometric parameters of the shell are taken to be R − = 0.289 m, L = 0.3048 m and h = 0.004 m or h = 0.04 m.
Due to symmetry of the problem, only one octant of the shell is modeled by fine regular meshes of 16 × 64 and 32 × 128 respectively for thick and thin cylindrical shells to describe correctly the through-thickness distributions of transverse stresses. Figures 12 and 13 present distributions of the transverse displacement of the bottom surface (z = −0.5) in section AB (θ 1 = 0) and the electric potential of the middle surface (z = 0) and the surface located close to the middle surface with z = −0.0555 in the same section, where z = θ 3 /h is the dimensionless thickness co-ordinate. A comparison with the isoparametric layer-wise piezoelectric shell element with linear interpolations of the displacements and electric potential in the thickness direction [50] is also given.
Tables 6 and 7 list the results of the convergence study due to increasing the number of SaS for thick and thin  piezoceramic cylindrical shells. Figures 14 and 15 show distributions of the transverse displacement, stresses, electric potential and electric displacement (48) through the thickness of the shell for different values of the material gradient index. These results demonstrate again the high potential of the GeXSaS4 element because the boundary conditions on bottom and top surfaces for the transverse stresses are satisfied properly. Figure 16 displays the results of the convergence study for a thin shell owing to mesh refinement through the normalized transverse displacement, electric potential, stresses and electric displacement for different material gradient indices. The reference values have been obtained taking five SaS and a regular 64 × 128 mesh. As can be seen, the GeXSaS4 element behaves well for thin shells in the case of coarse mesh configurations.   Figure 12: Distribution of the transverse displacement of the bottom surface u 3 (0, θ 2 , −0.5) and the electric potential of the middle surface ϕ(0, θ 2 , 0) in section AB for a thin cylindrical shell along coordinate θ 2 ; GeXSaS4 element (-) using five SaS and isoparametric layer-wise shell element [50] for α = 0 (∘). Figure 13: Distribution of the transverse displacement of the bottom surface u 3 (0, θ 2 , −0.5) and the electric potential of the surface located close to the middle surface ϕ(0, θ 2 , −0.0555) in section AB for a thick cylindrical shell along coordinate θ 2 ; GeXSaS4 element (-) using nine SaS and isoparametric layer-wise shell element [50] for α = 0 (∘).

Figure 14:
Through-thickness distributions of the transverse displacement, stresses, electric potential and electric displacement for a thick cylindrical shell using nine SaS.

Figure 15:
Through-thickness distributions of the transverse displacement, stresses, electric potential and electric displacement for a thin cylindrical shell using five SaS.

Figure 16:
Convergence study due to mesh refinement for a thin cylindrical shell using five SaS and regular 2n × 8n meshes with n = 1, 2, 4, 8, 12 and 16; reference solution is provided by a regular 64 × 128 mesh.

Conclusions
The paper presents a GeX hybrid-mixed FGM piezoelectric four-node solid-shell element based on the SaS formulation in which displacements and electric potentials of SaS are utilized as fundamental shell unknowns. The SaS are located at Chebyshev polynomial nodes inside the shell that improves the behavior of the higher-order Lagrange interpolations. To implement the efficient analytical integration throughout the element, the extended ANS method is employed to interpolate all components of the strain tensor and electric field vector. The element stiffness matrix of the proposed GeX piezoelectric solid-shell element has a correct rank and is evaluated without using the expensive numerical matrix inversion. The developed solid-shell element exhibits an excellent performance in the case of coarse meshes and can be recommended for the 3D coupled electro-mechanical analysis of thick and thin FGM piezoelectric doubly-curved shells.