Generalised assumed strain curved shell finite elements (CSFE-sh) with shifted-Lagrange and applications on N-T’s shells theory

Abstract We present a simple methodology to design curved shell finite elements based on Nzengwa-Tagne’s shell equations. The element has three degrees of freedom at each node. The displacements field of the element satisfies the exact requirement of rigid body modes in a ‘shifted-Lagrange’ polynomial basis. The element is based on independent strain assumption insofar as it is allowed by the compatibility equations. The element developed herein is first validated on analysis of benchmark problems involving a standard shell with simply supported edges. Examples illustrating the accuracy improvement are included in the analysis. It showed that reasonably accurate results were obtained even when using fewer elements compared to other shell elements. The element is then used to analyse spherical roof structures. The distribution of the various components of deflection is obtained. Furthermore, the effect of introducing concentrated load on a cylindrical clamped ends structure is investigated. It is found that the CSFE3-sh element considered is a very good candidate for the analysis of general shell structures in engineering practice in which the ratio h/R ranges between 1/1000 and 2/5.


Nomenclature a••, b••,
First, second and third fundamental c•• form of the surface S (aα , a 3) Local vector basis at the midsurface S a i (i = 1, ......, 9) Unknown coefficients a αβ , a αβ Covariant and contravariant components of the metric tensor of the surface S Aα Norm of aα b αβ Covariant component of curvature tensor b αβ b α β Contravariant and mixt components of curvature tensor Bm , B k , Bg Global matrix shape function of the change of first, second and third fundamental forms respectively (CB) ij = C ijkl B kl : Shifted-Lagrange space of shape functions Q•• Covariant linearized change of the third fundamental form tensor of S, or change of Gauss curvature tensor of the mid surface S uα , u 3 Components of a displacement in (aα , a 3) basis Uα , U 3 Components of the displacements in (gα , g 3) basiŝ︀ U Global discrete displacements of vector (ξα , ξ 3) , are displacement components (ηα , η 3) of vector fields on S

Introduction
A great number of research works has been expended over the past five decades on the development of shell finite elements methods for the analysis of curved structures [1,2]. The formulation of simple and robust finite elements has become one of the most important research fields in structural mechanics [3], like high order elements of Carrera et al. [4]; Keshava Kumar S. et al. [5]. However, surface displacement-based elements such as the four-node quadrilateral element behave very poorly for the case of bending problems [3]. Considerable efforts have been oriented to overcome the weaknesses of these elements by the development of efficient elements using different concepts and formulations such as the assumed strain or enhanced assumed strain elements Li and Huang [6]; Piltner and Taylor [7]; the generalized conforming elements Chen et al. [8]; Li and Huang [6]; the quasi-conforming elements Wang et al. [9]; Xia et al. [10] and the quadrilateral area coordinate elements Li and Huang [6]; Cen et al. [11]. Other robust membrane elements with in-plane rotation Kugler et al. [12]; Cen et al. [13]; Zouari et al. [14] have been developed. Grafton and Strome [15] developed conical segments for the analysis of shell of revolution. Later Jones and Strome [16] modified the method and used meridional elements which were found to lead to considerably improved results for the stresses. Curved rectangular and cylindrical shell elements were also developed by [17,18]. However, to model a spherical shape shells using the finite element method, triangular and rectangular spherical shell elements are needed [19][20][21]. Several spherical shell elements have been developed to analyse curved structures. Most notably, the higher order elements of Bathe et al. [21][22][23][24], Do-Nyun Kim et al. [2], Tornabene et al. [25], Zeighampour [26], Feumo et al. [27] and Echter et al. [28] have resulted in an improvement of the accuracy of the results. However, this improvement is achieved at the expense of increased computational efforts and storage. Meanwhile, a simple alternative strain based approach by István M. et al. [18], Lakhdar B. [29], Nkongho et al. [30] and Kim et al. [22] has been used to develop curved elements.
In this approach, the exact terms representing all the rigid body modes and the displacement functions representing the element's strain are determined by assuming independent strain functions insofar as it is allowed by the compatibility equations (A.I. Mousa et al. [19]). These elements were found to yield faster convergence when compared with other available finite elements Lakhdar B. [29], Nkongho et al. [30], Sabir et al. [31] and Mousa et al. [19,32].
The strain-based approach was also used to develop rectangular and triangular spherical shell elements. These spherical elements possess only the five essential external nodal degrees of freedom (dof) at each summit node and were found to have excellent convergence. Most spherical shell structures are supported by circular arched beams. In finite element analysis, these beams are usually modelled using finite elements having six degrees of freedom at each summit node to represent their general three dimensional forms of deformation. These six degrees of freedom are the usual five essential external degrees used in shell analysis together with a sixth representing the rotation about the normal to the shell surface Nkongho et al. [30] and István M. et al. [18]. It is therefore of interest to introduce in the shell an additional degree of freedom (as a sixth degree) representing this in plane (as drilling) rotation in order to enhance the compatibility between the spherical shell and curved beam elements Nkongho et al. [30]. Unlike, the strain based approach, also known as the Cardiff Approach (Mousa et al. [19,32]), in the present study, no rotational degree of freedom will be used on the N-T's shell equations [33].
The objective of this paper is to develop CSFE triangular shell finite elements that can be practically used for N-T shell equations which are more general for shell structures. In the following sections of the paper, we first present the formulation of general curved finite element in shifted-Lagrange space of shape functions. We then propose a simple methodology to design these elements using the strain approximation technique and apply it on an example. A large number of elements can be constructed using our methodology. We introduce a selected CSFE3-sh triangular shell finite element with its numerical test results.

Formulation of general curved
shell finite elements with shifted-La Grange (CSFE-sh)

N-T's two dimensional elastic thick shells model
The N-T's two dimensional model for linear elastic thick shells has been deduced from the three dimensional problem without any ad hoc assumption whether of a geometrical or mechanical nature.
Let M be a point in the domain Ω (see Figures 1 and 2), from its position we define the basis For a three dimensional shell, the following basis (︀ g α , g 3 )︀ ,(gα , g 3) are dual bases and (︀ a α , a 3 )︀ , (aα , a 3) defined on the midsurface also constitute dual bases of shell ifχ = h 2R < 1. h denotes the shell's thickness while its components defined as follow in the local coordinates system: (2)  We deduce from equations (1a) that g α = (︀ µ −1 )︀ α ν a ν and g 3 = a 3 .
The strain tensor for a three dimensional shell reads: From the limit analysis Nzengwa et al. [33] demonstrated that the displacements field satisfies the equation ϵ i3 (U) = 0 and the unique solution is (ξ α ) are local membrane displacement components and u 3 = ξ 3 the transvers displacement, The general shell deformation tensor derived from the kinematics defined above is: Where: e αβ : is the membrane strain tensor; k αβ : is the Bending strain tensor; Q αβ : is the Tensor of Gauss curvature (scarcely used in literature).

CSFE-sh shape functions space
The methodology is that of the strain-based approach.
Here the general strain tensor components are not simplified with respect to a particular shell geometry. So the interpolation here depends on general Assumed Strain using 'Shifted-Lagrange (SL)' shape functions basis. The strain tensor components of N-T shell model are function of membrane displacements ξα (x, y), their first order derivatives ξ α,β , transverse displacements ξ 3 (x, y), and its first and second order derivatives ξ 3,α and ξ 3,αβ R. Nzengwa [33]. For this reason, the strain tensor components can be approximated at least in IL p+2,q+2 because IL 1,1 functions cannot materialize the curvature and as consequence, final displacements results will certainly be compromised. Three independent unknown displacements are sufficient to express strain and stress tensors; then 3×n -coefficients are needed to totally define the field of displacements (n is the number of nodes). The two membrane displacements are in IL p+1,q+1 while the transverse displacement is in IL p+2,q+2 .

CSFE-sh model
As previously mentioned, the CSFE-sh model is built from strain assumption hypothesis on the N-T kinematic equations. The interpolation is based onC 0 functions, in-plane linear and transversally quadratic which must properly reproduce the rigid-body like motion. The Figure 4 below presents the shape of both CSFE3-sh and CSFE4-sh. In what follows, the curvilinear coordinates is replaced by (x, y, z). Let us recall that the two-dimensional displacements of the N-T shell kinematic equations reads and their consequences on the plane state deformations tensor reveal Here, e αβ , k αβ , Q αβ are respectively the membrane deformation tensor, the bending tensor and the Gaussian curvature tensor components which is not found in many shell theories [33].
The displacements which describe well the rigid-body like motion hold for ε αβ This can be verified if and only if e αβ (u0) = 0, k αβ (u0) = 0 and Q αβ (u0) = 0 (11) In what follows, the Einstein convention on repeated index is not applied to below equations. The strain tensor components defined in equation (11) yields Solutions to above partial differential equations (p.d.e) may be many but the simplest is Next, because of the 3 unknown displacements per node, we need 3 times n-independent coefficients to totally define the field of displacements over a given element. The final solution of the p.d.e issued from assumed strain can then be found. It means that, there exist functions of bilinear polynomials To make it simple, let us construct the functions fe (x, y) , f k (x, y) and f Q (x, y) using Shifted-Lagrange base shape functions. We therefore write Let u = (ξα (x, y) , ξ 3 (x, y)) be the final solution that satisfies equations (20), (21) and (22), there exist where IL 1,1 denotes the linear Shifted-Lagrange functions space and IL 2,2 is the quadratic Shifted-Lagrange functions space defined as follow: let be real coefficients, then X = αs x + βs and Y = αs y + βs are first order polynomial functions. IL p,q space corresponds to the set of monomers X i Y j with i = 0, 1, . . . , p and j = 0, 1, . . . , q. By utilizing Pascal monomers triangle for hierarchical triangular element or quadrangular element, this polynomial space can be presented as below.
The particular solution u d (ξ d α (x, y) , ξ d 3 (x, y)) previously mentioned is then presented according to the element type and such that the bilinear polynomial function G (fe (x, y) , f k (x, y) , f Q (x, y)) in equation (19) satisfies the following compatibility equations 1 2 The equation (23b) is obtained by considering the above nine components of strain [e11 e 22 e 12 k 11 k 22 k 12 Q 11 Q 22 Q 12] to be independent. As they are a function of the three displacements ξ 1 , ξ 2 and ξ 3 , they must satisfy three additional compatibility equations. These compatibility equations are derived by eliminating ξ 1 , ξ 2 and ξ 3 from Equations (10d), (10e) and (10f).

CSFE3-sh triangular element
CSFE3-sh is a double curved triangular shell finite element having three nodes and three unknown displacements per node; two membrane displacements ξ d α (x, y) and one transverse displacement ξ d 3 (x, y), which is different from that of Mousa [19] and Bathe [2,23] who included an additional rotational degree of freedom. The assumed strains components are approximated by polynomials that satisfy compatibility equations (23b) and resultant nodal displacements which are bilinear and quadratic shape functions in the space of shifted-Lagrange polynomials.
Let (x1, y 1) , (x2, y 2) and (x3, y 3) be curvilinear coordinates of nodes 1 , 2 and 3 respectively, a 1 , a 2 , a 3 , . . . , a 9 are the 3 × n -independent coefficients necessary to define the 3 × n -degree of freedom (see Figure 6) as below: ξ d 1 , ξ d 2 and ξ d 3 are displacement components which result from assumed strain functions with no energy contribution from rigid body like motion. Having used zero constant to capture the rigid body modes, the remaining nine constants are available for building the displacements field due to the strains within the element. These constants can be apportioned among the strains in several ways. For the present element, the following is proposed such that the assumed strains in a double curved triangular element CSFE3-sh hold for: Ax Ay Rx (a4 + a 5 Y + a 6 X) Ay Ry + 2Ay,x Ax Ay Ry (a1 + a 2 X + a 3 Y) The solutions fe (x, y) , f k (x, y) , f Q (x, y) satisfy very well the compatibility equations (23b). The resulting displacements in (24), (25) and (26) are shifted-Lagrange polynomials. The element has three degrees of freedom at each summit node and the rigid modes are exactly represented.

Variational problem
In Order to calculate the stiffness matrix, we have to formulate the variational problem over a domain. Let the border of S, ∂S = 0 ∪ 1 be partitioned in two parts and the border of the shell ∂Ω

}︁
. We suppose here that, the shell is clamped on Γo and loaded by volume and surface forces as stated above, the three dimensional variational equation related to the equilibrium reads: Where (Ω) ; ∂ j U i ∈ L 2 (Ω) and U i = 0 on Γ 0 } is the Sobolev space. : and . denote respectively tensors and vectors scalar products. p i = Volume forces in the domain Ω. p i = Surface forces on Γ.
Knowing that, the constitutive law of the linear elastic homogenous material is and C ijkl = λg ij g kl + µ (︁ g ik g jl + g il g jk )︁ , with λ = 2µλ λ + 2µ ; The problem is now stated as follow: The strain tensor ϵ = [︀ ϵ ij ]︀ is defined at (10c), if substitution is done in the left hand side of the equation (34), the problem is now equivalent to: The existence and uniqueness of the solution of the variational equations are found in Nzengwa et al. [33].
Naturally at the mid surface, g αβ = (︀ µ −1 )︀ α ρ (︀ µ −1 )︀ β λ a ρλ depends on z-parameter. By using Taylor's expansion on By truncating at n = 1 and integrating through the thickness we obtain the best first order two-dimensional equation: Where E and v are respectively the Young modulus and Poisson's coefficient for a given homogeneous elastic material. One can deduct from equation (37) what follows The right-hand side of the best first order twodimensional variational equation reads Here p i and p i are defined as in R. Nzengwa et al. [33]. Let F be the resultant force vector from distributed load F V + M such that: By equating A 1 (u, v) to L(v), we obtain the structural discrete equilibrium equation over the whole area as follows: Over a single curved element of area Se the elementary stiffness matrix Ke and elementary force vectors Fe are locally calculated as follow: )︂ dSe; P i e denotes element nodal concentrated load. By regular assembling process of elementary stiffness matrix Ke and elementary force vectors Fe, we then express the global stiffness matrix and force vector over the whole area as follow: Ke; (47) Fe .
Where N S is the total number of elements of the discretized domain.
A FEM model CSFE3-sh is built from N-T shell theory. The assumed strains are chosen in order to satisfy compatibility equations (23b) and resultant displacements are shape functions of shifted-Lagrange space of polynomials. The rigid body like motion of no independent coefficient is well described with zero energy contribution. Therefore, nine parameters are used to describe completely the behavior of a three nodes element. For the calculation of the stiffness matrix and the force vector related to the distributed load we formulate the total potential energy of a single element. We note that the expression above contains only the strain energy and the work of the distributed load. It is important to note that in this case, the term related to the concentrated loads is excluded, consequently the vector of concentrated loads P i e should be produced additionally. This is an easy task based on the nodal degrees of freedom. The advantage of the double curved shell finite element is that the curved surface is captured exactly; as a consequence, it provides fairly accurate results even if the number of elements is relatively small. The element developed herein is first tested by applying it to the solution of a benchmark shell problem [20,29], and is then used to analyse a clamped ends cylindrical structure. The distribution of the various components of deflection will be investigated.

Hemisphere under diametrically opposed loads
The test of a thin hemisphere (R/h = 250) subjected to the free base with four concentrated loadings (see Figure 7) is used to check the absence of membrane locking and the good representation of the rigid body like motion. In this example, the hemisphere undergoes significant rotations of rigid bodies around the normal of the mid surface. The inextensible bending deformations out of membrane are also significant and this problem serves consequently as an excellent test to examine the aptitude of a shell element to represent the rigid modes and the inextensible modes. A reference solution presented in [34] gives for displacement according to the direction of the load: U A = V B = 0.094mm. Because of symmetry considerations, only the quarter of the hemisphere is discretized. The mesh is regular and the number of elements varies from 2 to 20. We compare our results obtained with both CSFE3 and CSFE3-sh, to several finite elements of shells. In the Figure 8, the convergence of displacement U A according to the number of elements is represented, Table 1 presents the values of displacements U A obtained. It appears that the model CSFEsh with generalized Gaussian integration in the membrane and bending does not present any membrane locking but  Figure 9: Convergence of U A compared to other models very good convergence. So it is needless to decouple any of the behavior to improve neither our results nor using reduced integration.

Benchmark 2: Pinched cylinder
This test case makes it possible to check the behavior of the element in bending and shearing. This test was studied by several authors. The two ends of the shell are clamped by an infinitely rigid diaphragm. One quarter of the cylinder is meshed thanks to symmetries of the problem. The shell is subjected to a specific effort at point A, ‖P‖ = 1 ( Figure 10). The data of geometry, loading and of material for this test case are listed in Table 2.
This traditional test of cylinder subjected to two opposite loadings diametrically concentrated whose ends rest on two rigid diaphragms in their plan is of thin shell. It constitutes a severe case to study the capacity of a shell element to describe fields of deformations of complex membrane with a significant share of bending without extension of mid surface, in particular on the level of the zones where the forces are exerted. A reference solution was presented in Lindberg G.M [35]:    Displacement V D following Y:

Benchmark 3: Thick walled cylinder with free ends
One of the problems frequently used to evaluate the performances of finite element on thick shell theory is found in G.Raju [1]; A steel cylinder with homogeneous, isotropic Geometrically the entire cylinder is uniform (across the cross section also), material is isotropic in nature. Entire analysis work has been done assuming neglecting thermal effects. Due to symmetry of loads and geometry of the vessel, we model one quarter of the structure. All calculations and analysis are done at the mid surface of the thick cylinder.

Discussion
As compared to CSFE3 Nkongho et al. [30], which is a single curved triangular fem, the CSFE3-sh rate of conver-gence is high (see Figures 11 and 12). Because of its single curvature, CSFE3 cannot capture the geometry of a double curved structure that iss why it was not used to compute the hemisphere benchmark ( Figure 8 and Table 1). Several authors have solved the problem of hemisphere by using both triangular and rectangular finite elements, where all the displacement components are represented by cubic polynomials and also proposed a series (closed form) solution. Later on, Dawe used a quintet order triangular element, having 54 degrees of freedom to analyze the same shell as A.I. Mousa [19]. He reported that the results obtained were more accurate than those given by Yang for the same number of elements. However, this element has an excessively large number of degrees of freedom A.I. Mousa [19]. It would therefore require large computational effort to perform the analysis using this element. In comparison, the present element includes only three degrees of freedom.
The hemispherical shell described above is analyzed here using: the triangular spherical shell element that has only three degrees of freedom at each summit node. The results from the analyses are compared with those obtained from the series (closed form) solution given in A.I. Mousa [19,32] and Yang [36]. Convergence tests were carried out for the normal deflection at the point A (Figure 7) of the shell. Figure 9 shows that CSFE3-sh requires only 10×10 mesh size to converge to acceptable results with a difference of less than 3% in the case of the triangular element developed by Sabir and Djoudi [31,37,38] and less that 2.9% in the case of present element compared with the series solution, while it was reported that SFE3 element gives an error of above 10% for the same mesh size Mousa [19].
Further investigations on the deflection are shown in Figures 11, 12 and 13, which indicate excellent agreement between the results obtained from the present element and those found in the literature for the variation of normal deflecting along structure generic line. The pinched cylinder is a very severe benchmark to test a fem model (especially for thin shells). The loading case here is not symmetric (because of the concentrated point load) and the solution can't be handled by 3-D elasticity theory. The CSFE3sh shows quite good convergence for this membrane dominated shell problem. The convergence curves when the edges of the structure are clamped are shown. This is a difficult problem to solve when the thickness is small Mousa [19], but the problem is an excellent test case because of the negative Gaussian curvature of the shell mid surface. The element shows in fact good accuracy characteristics for the practical range of h/R=1/250; 1/100 and 2/5.
The last range h/R=2/5 falls under thick shells Nkongho et al. [30] (hollow shafts, spacers. . . ) and the accuracy is excellent than that of CSFE3. For only 6X6 mesh size, it presents an acceptable convergence results with a difference of less than 1.2% (Table 5) compared with the analytical solution from three dimensional elasticity theory. The CSFE3-sh appears to be a very good candidate to handle both thin shells and thick shells problems in engineering.

Conclusions
In this paper, we proposed a systematic procedure to construct spatially CSFE-sh triangular shell finite elements. The method is mechanically clear as well as simple and effective. We then constructed 3-nodes (CSFE3-sh) shell finite elements. For the selected element (the CSFE3-sh), we performed well-chosen numerical tests and showed convergence curves. While the elements have been developed and tested using the continuum-mechanics based approach with the N-T kinematics (with the underlying basic shell model identified by Nkongho et al. [30]), the same interpolation approach is of course also applied by Moussa [19] and Lakhdar [29]. The triangular family of elements considered is a good candidate for the analysis of general shell structures in engineering practice in which the range of h/R is usually between 1/1000 and 2/5. The elements show good behavior in the chosen test problems forth at range of thickness values. However, it is still necessary to study these elements further and to obtain uniformly optimal triangular shell finite elements that behave equally well in all types of shell problems.