Exact geometry solid-shell element based on a sampling surfaces technique for 3D stress analysis of doubly-curved composite shells

Abstract A hybrid-mixed ANS four-node shell element by using the sampling surfaces (SaS) technique is developed. The SaS formulation is based on choosing inside the nth layer In not equally spaced SaS parallel to the middle surface of the shell in order to introduce the displacements of these surfaces as basic shell variables. Such choice of unknowns with the consequent use of Lagrange polynomials of degree In − 1 in the thickness direction for each layer permits the presentation of the layered shell formulation in a very compact form. The SaS are located inside each layer at Chebyshev polynomial nodes that allows one to minimize uniformly the error due to the Lagrange interpolation. To implement the efficient analytical integration throughout the element, the enhanced ANS method is employed. The proposed hybrid-mixed four-node shell element is based on the Hu-Washizu variational equation and exhibits a superior performance in the case of coarse meshes. It could be useful for the 3D stress analysis of thick and thin doubly-curved shells since the SaS formulation gives the possibility to obtain numerical solutions with a prescribed accuracy, which asymptotically approach the exact solutions of elasticity as the number of SaS tends to infinity.


Introduction
A conventional way of developing the higher-order shell formulation consists in the expansion of displacements into power series with respect to the transverse coordinate, which is referred to the direction normal to the middle surface. For the approximate presentation of the displacement field, it is possible to use finite segments of power series because the principal purpose of the shell theory consists in the derivation of approximate solutions of elasticity. The idea of this approach can be traced back to Cauchy [1] and Poisson [2]. However, its implementation for thick shells is not straightforward since it is necessary to retain the large number of terms in corresponding expansions to obtain the comprehensive results.
An alternative way of developing the shell theory is to choose inside each layer a set of not equally spaced sampling surfaces (SaS) Ω (n) 1 , Ω (n)2 ,. . ., Ω (n)In parallel to the middle surface in order to introduce the displacement vectors u (n)1 , u (n)2 , u (n)In of these surfaces as basic shell variables, where In is the total number of SaS of the nth layer (In ≥ 3) and n = 1,2,. . .,N, where N is the number of layers. Such choice of displacements with the consequent use of the Lagrange polynomials of degree In − 1 in the thickness direction for each layer allows the presentation of governing equations of the layered shell formulation in a very compact form. The SaS concept was proposed recently by the authors to evaluate analytically the 3D stress state in rectangular plates and cylindrical and spherical shells [3][4][5].
It should be noticed that the SaS shell formulation with equally spaced SaS does not work properly with the Lagrange polynomials of high degree because of the Runge's phenomenon [6]. This phenomenon can yield the wild oscillation at the edges of the interval when the user deals with any specific functions. If the number of equispaced nodes is increased then the oscillations become even larger. However, the use of the Chebyshev polynomial nodes [7] inside the shell body can help to improve significantly the behavior of the Lagrange polynomials of high degree because such a choice permits one to minimize uniformly the error due to the Lagrange interpolation. This fact gives an opportunity to derive displacements and stresses with a prescribed accuracy employing the sufficiently large number of SaS. It means in turn that the solutions based on the SaS technique asymptotically approach the 3D exact solutions of elasticity as the number of SaS goes to infinity. Note also that the origins of the SaS concept can be found in contributions [8,9] in which three, four and five equally spaced SaS are employed. The SaS formulation with the arbitrary number of equispaced SaS is considered in [10]. The more general approach with the SaS located at the Chebyshev polynomial nodes was developed later [3][4][5]. In contribution [11], both SaS formulations are compared with each other.
It is well-known that the low-order finite elements are sensitive to shear locking. This is because of incorrect shear modes, which infect the pure bending element behavior. In regards to exact geometry or geometrically exact (GeX) four-node shell elements, they reproduce additionally membrane locking. The abbreviation "GeX" reflects the fact that the parametrization of the middle surface is known a priori and, therefore, the coefficients of the first and second fundamental forms are taken exactly at element nodes [12,13]. The feature of GeX solid-shell elements [12,13] is that they are based on strain-displacement relationships of the simplest SaS formulation with In = 3, which precisely represent all rigid-body motions of the shell in any convected curvilinear coordinates. This fact is of great importance since one may read in paper [14] that "shell theory is an absolute academic exercise" due to "the difficulties of representing the rigid body modes in shell finite element formulations".
To avoid shear and membrane locking and have no spurious zero energy modes, the hybrid-mixed finite element method can be applied. This method pioneered by Pian [15] utilizes the robust interpolation of displacements on the element boundary to provide displacement compatibility between elements, whereas the internal stresses are assumed so as to satisfy the differential equilibrium equations. The Pian's hybrid stress finite element was originally based upon the principle of the stationary complementary energy. Alternatively, the assumed stress finite element was proposed by applying the Hellinger-Reissner variational principle that simplifies the evaluation of the element stiffness matrix [16]. Then, the assumed strain [17] and assumed stress-strain [18] finite elements were developed. The former is based on the modified Hellinger-Reissner functional in which displacements and strains are utilized as basic shell variables, whereas the latter departs from the Hu-Washizu functional depending on dis-placements, strains and stresses. However, here we do not use these terms because a more general term hybrid-mixed finite element covers all hybrid and mixed finite elements, which are "formulated by multivariable variational functional, yet the resulting matrix equations consist of only the nodal values of displacements as unknowns" [19].
In the present paper, the hybrid-mixed GeX four-node solid-shell element formulation is proposed in which the SaS are located inside each layer at Chebyshev polynomial nodes [5]. To circumvent locking phenomena, the assumed interpolations of displacement-independent strains and stress resultants are invoked and utilized together with displacement and displacement-dependent strain interpolations into the Hu-Washizu variational equation. Such an approach exhibits an excellent performance in the case of coarse mesh configurations and has computational advantages compared to conventional isoparametric hybridmixed solid-shell element formulations [20][21][22][23] because it reduces the computational cost of numerical integration in the evaluation of the element stiffness matrix. This is due to the facts 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 shell element formulation. Second, the GeX four-node solid-shell element formulation is based on the effective analytical integration throughout the finite element by the use of the enhanced ANS method [24,25]. The latter has a great meaning for the numerical modelling of thick doubly-curved shells with variable curvatures.
The thin doubly-curved shells based on strong and weak formulations are widely discussed in the literature. The static and dynamic analyses of doubly-curved shells by using the classic Kirchhoff-Love and Timoshenko-Mindlin-type theories can be found in books [26][27][28][29][30]. The state-of-the-art development of the problem is analyzed in review articles [31,32]. The higher-order doubly-curved shell formulation based on the Carrera's equivalent single layer theory [33] accounting for thickness stretching has been proposed in [34,35]. Both free vibration and static problems are discussed with a particular emphasis on the stress recovery procedure. The authors report that their procedure leads to stable, accurate and reliable results for the moderately thick and thin doubly-curved shells with variable principal curvatures. However, for the analysis of thick doubly-curved shell structures instead of the postprocessing stress recovery technique a more general approach based on the 3D constitutive equations should be applied. Such a question is discussed here in detail. ribe all SaS of the nth layer and run from 1 to n I ; Latin tensorial indices l k j i , , , e from 1 to 3; Greek indices  , range from 1 to 2.   [7]. This fact has a great meaning for a ergence of the SaS method [3,11].
inematic description of deformed shell sition vector of the deformed shell is written as

Kinematic description of undeformed shell
Consider a layered 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; r = r(θ 1 , θ 2 ) is the position vector of any point of the middle surface; a i (θ 1 , θ 2 ) are the base vectors of the middle surface given by where c (n)in are the components of the shifter tensor at SaS.
Here and in the following developments, (. . .) ,i stands for the partial derivatives with respect to coordinates θ i ; the index mn identifies the belonging of any quantity to the inner SaS of the nth layer and runs from 2 to In − 1, whereas the indices in, jn, kn describe all SaS of the nth layer and run from 1 to In; Latin tensorial indices i, j, k, l range from 1 to 3; Greek indices α, β range from 1 to 2.
Remark 1: As follows from (3), the transverse coordinates of inner SaS coincide with the nodes of Chebyshev polynomials [7]. This fact has a great meaning for a convergence of the SaS method [3,11].

Kinematic description of deformed shell
A position vector of the deformed shell is written as where u is the displacement vector, which is measured in accordance with the total Lagrangian formulation from the initial configuration to the current configuration directly.
In particular, the position vectors of SaS of the nth layer areR (n)in = R (n)in + u (n)in , where u (n)in (θ1, θ 2) are the displacement vectors of SaS of the nth layer Ω (n)in . The base vectors in the current shell configuration are defined asḡ In particular, the base vectors of deformed SaS of the nth layer (see Figure 2) arē where β (n)in (θ1, θ 2) are the values of the derivative of the displacement vector with respect to the thickness coordinate at SaS of the nth layer.  The Green-Lagrange strain tensor in an orthogonal curvilinear coordinate system [12] can be written as where A 3 = 1 and c 3 = 1. In particular, the Green-Lagrange strain tensor at SaS is where ε (n)in ij (θ1, θ 2) are the strains of SaS of the nth layer. Substituting (5) and (10) into strain-displacement relationships (13) and discarding the non-linear terms, one derives Next, we represent the displacement vectors u (n)in and β (n)in in the middle surface frame e i as follows: Using (15) and formulas for the derivatives of unit vectors e i with respect to orthogonal curvilinear coordinates [25], we have 1 Aα where Substitution of (16) and (17) into strain-displacement relationships (14) yields the component form of these relationships

Displacement and strain distributions in thickness direction
Up to this moment, no assumptions concerning displacement and strain fields have been made. We start now with the first fundamental assumption of the proposed higher order layer-wise shell theory. Let us assume that the displacements are distributed through the thickness of the nth layer as follows: where L (n)in (θ3) are the Lagrange polynomials of degree In − 1 expressed as The use of (11), (16) and (20) yields where M (n)jn = L (n)jn ,3 are the derivatives of Lagrange polynomials. The values of these derivatives at SaS are calculated as for jn ≠ in, Thus, the key functions β (n)in i of the proposed layer-wise shell theory are represented according to (22) as a linear combination of displacements of SaS of the nth layer u (n)jn i . Geometry solid-shell element based on analysis of composite shells

5
The following step consists in a choice of the consistent approximation of strains through the thickness of the nth layer. It is apparent that the strain distribution should be chosen similar to the displacement distribution (20), that is, The proof of this proposition follows from a homogeneous system of linear equations ...
and the identity ∑︁

Hu-Washizu mixed variational equation
To develop the assumed stress-strain finite element formulation, we invoke the Hu-Washizu variational principle in which displacements, strains and stresses are utilized as independent variables [36]. It can be written as follows: where σ (n) ij are the stresses of the nth layer; e (n) ij are the displacement-independent strains of the nth layer; C (n) ijkl are the elastic constants of the nth layer; u − i and u + i are the displacements of bottom and top surfaces; p − i and p + i are the tractions acting on the bottom and top surfaces; c − α = 1− kαh/2 and c + α = 1+ kαh/2 are components of the shifter tensor of bottom and top surfaces; W Σ is the work done by external loads applied to the edge surface Σ. Here and in the following developments, the summation on repeated Latin indices is implied.
Following the SaS technique, we introduce the third assumption of the proposed hybrid-stress solid-shell element formulation. Let the displacement-independent strains be distributed through the thickness similar to displacement and displacement-dependent strain distributions (20) and (24), that is, where e (n)in ij are the displacement-independent strains of SaS of the nth layer. Substituting strain distributions (24) and (30) in (28) and (29), and introducing stress resultants the following variational equation is obtained:

Hybrid-mixed GeX solid-shell element formulation
The finite element formulation is based on the simple and efficient interpolation of shells via curved GeX four-node solid-shell elements  where Nr (ξ1, ξ 2) are the bilinear shape functions of the element; u (n)in ir are the displacements of SaS Ω (n)in at element nodes; ξα = (θα − dα) /ℓα are the normalized curvilinear coordinates (see Figure 3); 2ℓα are the lengths of the element in (θ1, θ 2) -space; the nodal index r runs from 1 to 4. The surface traction vector is also assumed to vary bilinearly throughout the element.
To implement analytical integration throughout the element, we employ the enhanced ANS interpolation [24,25] where ε (n)in ijr are the strains of SaS of the nth layer at element nodes. These strains can be evaluated as Here, B (n)in  0, 1, . . . , N) ,

Remark 3:
The main idea of such approach can be traced back to the ANS method [37][38][39] developed by many scientists for the isoparametric displacement-based and hybrid-mixed finite element formulations [22,23,28,40,41]. In contrast with above formulations, we treat the term "ANS" in a broader sense. In the proposed GeX four-node solid-shell element formulation, all components of the strain tensor are assumed to vary bilinearly inside the biunit square. This implies that instead of the expected nonlinear interpolation due to variable curvatures in straindisplacement relationships (18) and (19) the more suitable bilinear ANS interpolation is utilized.
Remark 4: In order to circumvent curvature thickness locking for the isoparametric four-node solid-shell element, Betsch and Stein [42] proposed to apply the bilinear interpolation (36) only for the transverse normal strain. It is apparent that curvature thickness locking is not related to the curved GeX four-node solid-shell element because it can handle the arbitrary geometry of surfaces properly. We advocate the use of the enhanced ANS method for all components of the strain tensor to implement the efficient analytical integration throughout the element.
From the computational point of view it is convenient to represent the ANS interpolation (36) as follows: where Here and below, the indices r 1 and r 2 run from 0 to 1. To overcome shear and membrane locking and have no spurious zero energy modes, the robust displacementindependent strain and stress resultant interpolations are utilized  where Qr 1 r2 are the projective matrices given by Inserting interpolations (34), (39), (42) and (43) into the Hu-Washizu variational equation (32) and replacing the metric product A 1 A 2 in surface integrals by its value at the element center, one can integrate analytically throughout the finite element. As a result, the following equilibrium equations of the hybrid-mixed GeX solid-shell element are obtained: e (n)in r1 r2 = (Qr 1 r2 ) T B (n)in r1 r2 U for r 1 + r 2 < 2, where F is the element-wise surface traction vector. Eliminating column matrices e (n)in r1 r2 and H (n)in r1 r2 from (45)-(47), we arrive at the governing finite element equations where K is the element stiffness matrix of order 12N SaS × 12N SaS defined as   Due to symmetry, only one quarter of the plate with dimensions 1  a and 3  b is discretized by uniform meshes of SaSGeX4 elements depicted in Figure 5. To compare the results with Pagano's exact solution [43], we introduce the following dimensionless variables at crucial points as functions of the dimensionless thickness coordinate z:  The data listed in Tables 1 and 2 show that the SaSGeX4 element allows reproducing the exact solution of elasticity [43] for thick and thin plates with a high accuracy by using a fine 64 64  mesh and the sufficiently large number of SaS. As can be seen, the SaSGeX4 element provides from four to five right digits for displacements and stresses except for the transverse normal stress comparing to authors' exact SaS solution [4]. Figure 6 displays the distributions of displacements and stresses (51) through the thickness of the plate for different values of the slenderness ratio h a / It is worth noting that the element stiffness matrix (49) requires only direct substitutions, i.e., no expensive matrix inversion is needed to derive it. This is impossible for the isoparametric hybrid-mixed finite element formulations [20][21][22][23]. Furthermore, the stiffness matrix is evaluated by using analytical integration throughout the element that permits the use of coarse mesh configurations. Thus, the hybrid-mixed GeX solid-shell element formulation developed is economical and efficient.

Numerical examples
The performance of the proposed hybrid-mixed GeX fournode solid-shell element on the basis of the SaS technique, denoted here by the SaSGeX4 element, is assessed through the use of exact solutions of elasticity for the laminated composite rectangular plate and spherical shell, and authors' hyperbolic composite shell example.

Laminated composite rectangular plate under sinusoidal loading
Here, we study a simply supported laminated composite rectangular plate subjected to a sinusoidally distributed transverse load acting on its top surface where a and b are the length and width of the plate (see Figure 4) Due to symmetry, only one quarter of the plate with dimensions a = 1 and b = 3 is discretized by uniform meshes of SaSGeX4 elements depicted in Figure 5. To compare the results with Pagano's exact solution [43], we introduce the following dimensionless variables at crucial points as functions of the dimensionless thickness coordinate z: The data listed in Tables 1 and 2 show that the SaS-GeX4 element allows reproducing the exact solution of elasticity [43] for thick and thin plates with a high accuracy by using a fine 64 × 64 mesh and the sufficiently large number of SaS. As can be seen, the SaSGeX4 element provides from four to five right digits for displacements and stresses except for the transverse normal stress comparing to authors' exact SaS solution [4]. Figure 6 displays the distributions of displacements and stresses (51) through the thickness of the plate for different values of the slenderness ratio a/h utilizing five SaS inside each layer and the  same 64 × 64 mesh. These results demonstrate convincingly the high potential of the SaSGeX4 element because the boundary conditions on bottom and top surfaces and the continuity conditions at layer interfaces for transverse stresses are satisfied correctly without employing the postprocessing stress recovery technique. Figure 7 shows the results of the convergence study due to mesh refinement through the use of normalized displacements and stresses for slenderness ratios a/h = 4 and 100 by choosing five SaS for each layer. The analytical answers are provided again by the exact SaS solution [4]. It is seen that the SaSGeX4 element behaves practically insensitive with respect to the mesh parameter n introduced in Figure 5.

Hyperbolic composite shell under nonuniform inner pressure
Finally, we study a hyperbolic composite shell subjected to inner pressure where ) , ( P 2 1   is the point belonging to the middle surface.
Due to symmetry of the problem, only one sixteenth of the shell is discretized with regular meshes depicted in Figure 10. Figure

Spherical shell under uniform inner pressure
Next, we consider a single-layer spherical shell with 0.02 ∘ hole at the top subjected to uniform inner pressure p − 3 = −p 0 . This problem is a good benchmark to test the proposed analytical integration scheme and verify that the SaSGeX4 solidshell element is free of locking. In the literature there is Lame's closed-form solution [44,45], which can be written as where r is the radial distance from a point to the origin; R is the radius of the middle surface. The mechanical and geometrical parameters of the shell are chosen to be E = 10 5 , ν = 0.3, R = 10 and ϑ * = 89.98 ∘ , where ϑ * is the up limit value of the meridional coordinate θ 1 (see Figure 8). To compare the results derived with Lame's analytical solution, we introduce dimensionless variables u 3 = 10Ehu 3 (P, z)/R 2 p 0 ,σ 11 = 10hσ 11 (P, z)/Rp 0 ,σ 33 = σ 33 (P, z)/p 0 , ε 11 = 10hEε 11 (P, z)/Rp 0 ,ε 33 = 10hEε 33 (P, z)/Rp 0 , z = θ 3 /h,  where P(0, 0) is the point belonging to the middle surface.
Owing to symmetry, we consider a part of the shell and model it by regular meshes shown in Figure 8. Tables 3 and 4 list the results for thick and thin spherical shells due to increasing the number of equispaced SaS by using a fine 64 × 1 mesh of SaSGeX4 elements. A comparison with Lamé's solution [45] is also given. As turned out, the SaSGeX4 element provides from 3 to 4 right digits for displacements and stresses. Figure 9 displays the through-thickness distributions of stresses for different values of the slenderness ratio R/h choosing seven SaS inside the shell body. One can see that the boundary conditions for the transverse normal stress on outer surfaces are satisfied properly. Table 5 lists the results of the convergence study due to mesh refinement for a thin shell with three equispaced SaS. It is seen that the SaSGeX4 element is free of locking and practically insensitive to coarse meshes.

Hyperbolic composite shell under nonuniform inner pressure
Finally, we study a hyperbolic composite shell subjected to inner pressure p − 3 = −p 0 cos 4θ 1 . The mechanical and geometrical parameters of the shell are taken as E 1 = 40E, E 2 = E 3 = E, G 12 = G 13 = G 23 = 0.6E, E = 10 6 , ν 12 = ν 13 = ν 23 = 0.25, r = 7.5, R = 15 and a = 20, where 2a is the height of the shell (see Figure 10). To analyze the numerical results efficiently, we introduce dimensionless stresses as functions of the dimensionless thickness coordinate z = θ 3 /h as follows: σαα = 10a 2 σαα(P, z)/h 2 p 0 , σ α3 = 10aσ α3 (P, z)/hp 0 , σ 33 = σ 33 (P, z)/p 0 , where P(θ 1 , θ 2 ) is the point belonging to the middle surface. Due to symmetry of the problem, only one sixteenth of the shell is discretized with regular meshes depicted in Figure 10. Figure 11 presents the through-thickness distributions of stresses (54) for different values of the slenderness ratio R/h by choosing seven SaS inside the shell body. One can see that the boundary conditions for transverse stresses on the bottom and top surfaces are satisfied again correctly without the use of the post-processing stress recovery technique.

Conclusions
The paper presents a geometrically exact hybrid-mixed four-node solid-shell element SaSGeX4 based on the SaS formulation in which the displacements of SaS are utilized as fundamental shell unknowns. The SaS are located at Chebyshev polynomial nodes inside the shell body that allows one to minimize uniformly the error due to Lagrange interpolations of displacements and strains through the thickness. To implement the efficient analytical integration throughout the element, the enhanced ANS method for all components of the strain tensor is employed. The element stiffness matrix is evaluated without the use of expensive numerical matrix inversion and has six zero eigenvalues as required for satisfaction of the general rigid-body motion representation. The SaSGeX4 element exhibits a superior performance in the case of coarse mesh configurations for all 3D benchmarks considered. It can be recommended for the 3D stress analysis of thick and thin doublycurved shells because the SaS solutions asymptotically approach 3D solutions as the number of SaS tends to infinity.