Finite Element Model for Linear Elastic Thick Shells Using Gradient Recovery Method

This research purposed a new family of finite elements for spherical thick shell based on Nzengwa-Tagne’s model proposed in 1999. The model referred to hereafter as N-T model contains the classical Kirchhoff-Love (K-L) kinematic with additional terms related to the third fundamental form governing strain energy. Transverse shear stresses are computed and C0 finite element is proposed for numerical implementation. However, using straight line triangular elements does not guarantee a correct computation of stress across common edges of adjacent elements because of gradient jumps.The gradient recoverymethod known as Polynomial Preserving Recovery (PPR) is used for local interpolation and applied on a hemisphere under diametrically opposite charges. A good agreement of convergence results is observed; numerical results are compared to other results obtained with the classical K-L thin shell theory. Moreover, simulation on increasing values of the ratio of the shell shows impact of the N-T model especially on transverse stresses because of the significant energy contribution due to the third fundamental form tensor present in the kinematics of this model. The analysis of the thickness ratio shows difference between the classical K-L theory and N-T model when the ratio is greater than 0.099.


Introduction
N-T's theoretical approach which was mathematically and rigorously deduced from three-dimensional linear elastic curvilinear media through multiple scaling and limit analysis is a more general Kirchhoff-Love (K-L) model.The displacement here is a two-degree polynomial of the thickness parameter  while the strain tensor which is planar contains the change in the third fundamental form in addition to the change of the first and second.
For more than thirty years, numerous articles, books, and theses have addressed the problem of shell.Plates and more generally thin shells represent over 70% of industrial calculations [1,2].The mechanical models have been validated by some well-known benchmarks.Some locally stiffened thin shells or more generally thick shells have not received the same development probably because some few existing models do not account for transverse stresses and have not been mathematically established.Moreover it is well known that transverse stresses can not be neglected as the shell becomes thicker.Without any ad hoc geometrical hypothesis, it was deduced that the strain tensor is planar; that is, ∈ 3 = 0 and, as a consequence, the limit displacement reads  =     +  3  3 , with   (, ) =   () −   () +  2   (),  3 =  3 () (where { 1 ,  2 ,  3 } is the three-dimensional contravariant shell basis); the in-plane strain ∈  =   −   +  2   and the stresses   depend on the in-plane strains and  3 are computed from appropriate differential equations; see [3].
The form of the displacement clearly shows that  0 finite element is discontinuous across adjacent elements and usually provides inaccurate results at elements boundaries.Some authors indeed proposed different numerical methods: the finite difference method (FDM); finite volume method (FVM); finite element method for various simulations like magnetohydrodynamic (MHD) [4][5][6]; Curved Triangular Finite Elements (CTFE) of [7]; a weighted average method to calculate stresses, which provided good results for both interior and boundary elements; and  2 -projection to also calculate stresses, by considering "discrete smoothing" and least squares fitting at the gauss points.But this last method was limited to one element; consequently, the smoothen stress is still discontinuous across element boundaries.This problem was completely solved 20 years ago, when Sheikholeslami et al. introduced their Superconvergent Patch Recovery (SPR) [8][9][10], where the discrete least square fitting was performed on an element patch, a set of elements having the same vertex.This SPR method produces a continuous stress field, which is superconvergent under uniform mesh.Soon after, Wiberg et al. incorporated equilibrium and boundary conditions to enhance SPR [11,12] and discussed strategies to improve the finite element solution  ℎ itself (other than the stress, which is essentially the gradient of  ℎ ) [12].Recently, Naga and Zhang and their colleagues proposed an alternative strategy, called Polynomial Preserving Recovery (PPR) [13][14][15][16][17], to recover the gradient.Theoretical analyses revealed that PPR has better superconvergence (over any mesh) properties than SPR and the numerical tests indicated that the a posteriori error estimator based on PPR is as good as or better than that of SPR ( [18,19], remarks: page 323).
In order to validate the N-T model, the program is first tested successfully on the widely known hemisphere under diametrically opposite charges of thin shells.Then thin shell theory of K-L and thick shell theory of N-T are implemented in order to evaluate the impact of thickness ratio.Next, simulations with various values of the characteristic shell parameter (thickness ratio) are implemented in order to reveal the contribution of Gauss curvature (change of the third fundamental form) in the stiffness energy.
In addition to that, Section 2 presents materials and methods: a brief description of the N-T model is presented without all the mathematical development, including gradient recovery method.A variational formulation of the shell equation and the resolution of transverse stresses equations are done.Section 3 is devoted to the discretization of the N-T model.Finite element spaces are next described and also the discretization scheme is layout.All numerical integrations are performed on a reference triangular element using the gradient recovery PPR method.Section 4 is completely devoted to the validation of the finite element model.Finally, results are discussed and concluded.

The N-T Model of Thick Shells
is the thickness, and  is the coordinate of  in ) denote a shell.We assume the surface  is bounded and sufficiently smooth for all subsequent computations.Let { 1 ,  2 ,  3 } and { 1 ,  2 ,  3 } denote the covariant and contravariant basis of the midsurface and { 1 ,  2 ,  3 } and { 1 ,  2 ,  3 }, respectively, the covariant and contravariant basis of the shell.Then where    =     and   denote curvature tensor components and   is the contravariant component of the metric of the midsurface .The repeated index convention is adopted.Values of Greek indices ,  take range in the set {1, 2} while Latin indices ,  take their values in the set {1, 2, 3}.A vector field can be expressed component wise indifferently in  basis or   -basis as follows: Then the strain tensor (see [3]) is given by (/ and ∇ indicate covariant derivation in Ω and , resp., while  , = /  ).The equation ∈ 3 () = 0 yields the following results: The strain tensor now reads Let us recall that (  ), (  ), and (  ) are, respectively, changes in the first, second, and third fundamental forms tensors, while the displacement is a more generalized Kirchhoff-Love displacement [3] and can also be written in Reissner-Mindlin format as follows: where   =  3, + 2 ] ] =   +  ]   ] and   are the rotations angles.Strain tensors are miscalculated in some literatures because wrong basis vectors were used.The condition ∈ 3 () = 0 modifies the constitutive law which, expressed with the Lamé constants, now reads (8) where  = 2/( + 2), or equivalently with Young's modulus  and Poisson's coefficient ], reads Consider a shell of thickness ℎ, clamped on a part of its border Γ 0 , subject to volume forces   and  3 and to surface forces ℎ  and ℎ 3 on the rest of its border Γ 1 .Suppose the forces are sufficiently smooth; then the transverse stresses are solutions to the differential equations: where   ∈  2 (Ω);  3 ∈  1 (−ℎ/2, ℎ/2;  −1 ()) and where  3 ∈  1 (Ω);  33 ∈  2 (−ℎ/2, ℎ/2;  −2 ()).
The generic point of the sphere is described by where  and  are curvilinear coordinates, 0 ≤  ≤ 2Π, and 0 ≤  ≤ Π;  is the radius of the sphere, and , ,  are the global coordinates.
The covariant and contravariant metric tensors on the middle surface  are defined by The covariant and contravariant curvature tensors on  (second fundamental form) are given by Let ℎ be the thickness of the sphere; the Christoffel symbols are defined by 2.2.Variational Equations.Let the border of ,  =  0 ∪ 1 be partitioned in two parts and the border of the shell Ω = Γ 0 ∪ Γ 1 with Γ 0 =  0 ×{−ℎ/2, ℎ/2}, and We denote Γ − =  × {−ℎ/2} and Γ + =  × {ℎ/2}.Suppose the shell is clamped on Γ 0 and subject to volume and surface forces as stated above; then the three-dimensional variational equation related to the equilibrium equation reads where and   = 0 on Γ 0 } is the Sobolev space and : and ⋅ denote, respectively, tensors and vectors scalar products.The variational formulation under Kirchhoff-Love's approach which is given by is a truncated thick shell or the best first-order thick shell variational equation under N-T's model; see [3].It can be observed that this equation has similarities with familiar equations in engineering literature.Let Nzengwa and Tagne Simo established existence and unicity of solutions of this truncated problem in  ad .We should remember that the displacement calculated after  is or in the Nzengwa-Tagne format which suggests  0 finite element implementation.
The linearized membrane strain tensor is expressed as follows: Using the formulas we can put (22) in the following form: where where (26) Mathematical Problems in Engineering 5 (, ) ;  2  (, 2) ; −    ; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0] ; The linearized part of the change of the curvature tensor of the middle surface is expressed as follows: The linearized part of the change of the third fundamental form tensor can be expressed as follows: (29)

Resolution of the Transverse Stress Ordinary Differential
Equation.Using expressions ( 27), (28), and (29), let  = −4/(ℎ − ), Γ   be the three dimension Christoffel symbols, and Then the first shear stress of ( 10) is given as follows: In the same way, let The second shear stress (10) is defined as follows: Let Then the third expression of ( 11) is defined as follows:

Polynomial Preserving Recovery Method.
Let us recall that the linear stress-strain relation does not guarantee smoothness across elements because of gradient jumps since the displacement is only  0 .The shell is meshed with straight lines triangular elements obtained by linear transformation of a reference 2D triangle.Jumps across elements will be prevented by implementing the gradient recovery methods that we briefly present hereafter.Let  ℎ be  0 finite element approximation of the solution , and let be two triangular elements with a common node .We denote are the vector of nodal values in  1 and  2 , respectively.Let   be the shape functions matrix; then restriction of  ℎ in each element reads where    stands for the derivative according to  or  of   .At the common node , the derivative of  ℎ ,   ℎ (), is not necessarily the same in both elements.This means that on a patch including a node  common to different elements the gradient of  ℎ , ∇ ℎ , cannot be calculated because of the jump on each part of the patch.Consequently ∇ ℎ is not a good approximation of ∇ across elements and stress and strain calculated using ∇ cannot be continuous across edges.How to define a unique ∇ ℎ at a node  common to different elements has been addressed by gradient recovery methods.The idea is to define a local operator  ℎ such that  ℎ  ℎ () is unique through any patch and | ℎ  ℎ − ∇| has better approximation than |∇ ℎ − ∇|.Consequently convergence issues of these methods should be considered also.In this analysis we consider the 2D PPR gradient recovery methods which also guarantees a superconvergence property of  ℎ  ℎ to ∇ independently of mesh size.The method is described as follows.
Let   be a node where ∇ is to be determined.Let    be the nodes of all triangles   having   as common vertex.Suppose  ℎ|  ∈  +1 , the set of polynomials of order  + 1; then the gradient recovery PPR operator consists in defining () ∈  +1 such that Mathematical Problems in Engineering 7 This value also depends on the sampling points chosen.In the 2D case it is proved (cf.[14]) that if  ≧  = ( + 2)( + 3)/2 and if the sum of two adjacent angles in the mesh is not more than , then  ℎ  ℎ (  ) is unique for any   .In this work the angle condition will be implemented in the triangularization.
The number  of sampling nodes will be respected as follows: For the 2D case, three types of nodes can be distinguished: internal nodes, boundary nodes at a corner, and boundary nodes out of a corner.
(i) For the in-plane displacement   ( = 1; 2), the existence and uniqueness solution are possible if ( 1 ,  2 ) ∈ ( 2 ()) 2 ; we fit linear polynômial using the same regular pattern: We scale by a factor ℎ with  = ℎ and  = ℎ, (, ) with respect to the six derivative values at the barycentric center of each element on the patch.Now we define  = ( Let    be the values of   at the nodal points     = 1, 2, 3, 4, 5, 6; we denote  = (   ) 6 =1 ; we hereby determine â so that â =  ℎ .
Then, we obtain notice that For the transversal displacement,  3 , the existence and uniqueness of a gradient recovery solution are possible if ( 3 ) ∈  2 (); we fit quadratic polynômial using the same regular pattern: ) and consider the diagonal matrix  1 = diag(1, ℎ, ℎ, ℎ 2 , ℎ 2 , ℎ 2 ) and also  ℎ = ( 30 ,  31 ,  32 ,  33 ,  34 ,  35 ,  36 ) the approximation vector of  3 in a nodal point   .Then we fit  2 (, ) = (1, , ,  2 , ,  2 ) −1 1 (  ) −1    ℎ and obtain the recovered gradient at the patch center point as follows: With  ℎ   given at each vertex by the same processes in (43) and (45), we are able to form a recovered gradient field by using the finite element basis functions.Recovering the gradient at a boundary vertex is more delicate.Using efficient strategy computational experiment indicated in [14], to recover the gradient at a vertex   ∈ , we look for the nearest layer of vertices around   that contain at least one internal vertex.Let this layer be the th one and denote the internal vertices in this layer by  1 ,  2 , . . .,   , where  ≥ 1.The union of the sampling points used in recovering the gradient at  1 ,  2 , . . .,   and the mesh nodes in the first  layers around   constitute the set of sampling point for recovering the gradient at   ; see [15].

Discretization of N-T Model
3.1.Finite Element Space.The continuous truncated variational equation ( 18) is defined in the space  ad .In order to define a finite element space, we begin by recalling the following results.Let   be a triangle with nodes  1 = ( 1 ,  1 ),  2 = ( 2 ,  2 ), and  3 = ( 3 ,  3 ).Let  1 and  2 be the sets of the first-and second-order polynomials with basis {1, , } and {1, , ,  2 ,  2 , }, respectively, then { 1 ,  2 ,  3 } and { 1 ,  2 ,  3 , 4 1  2 , 4 1  3 , 4 2  3 } also generate  1 and  2 where the barycentric basis functions of the triangle are defined by Let  ℎ be a triangularization of the midsurface of the shell and  ℎ the number of triangles.We denote by the subspace of dimension  ℎ (which is the number of nodes) and is the closure of , the subspace of dimension  2 ℎ which is the number of unisolvent points (nodes and edge mid points).
The space In the above formula,  ℎ is the gradient operator and  4 3 ,  5  3 , and  6 3 are, respectively, the unknown values of  3 in the midedges  1  2 ,  1  3 , and the coefficient of the approximation of the gradient when the vertex is internal of the mesh, the coefficient of the approximation of the gradient when the vertex is in the boundary at the corner of the mesh, the coefficient of the approximation of the gradient when the vertex is in the boundary not at the corner of the mesh.
Exploiting the above, let us take linear element on uniform triangular mesh of a regular pattern; we hereby investigate the case of the element where all three nodes are internal of the mesh; the matrices deduced from the discretization of the gradient are given as follows: When we return to energy function with  = 0, 1, we have Using ( 43), (45), and (51) in which all three nodes of the triangular element are internal, let  = 1, 2, 3; ,  = 1, 2 and we consider Then the deformation vector can be written as { Mathematical Problems in Engineering 3.3.Stiffness Matrix.Using local element, stiffness matrix can be expressed as follows: where is the generalized 19 × 19 behavior matrix.
Using (18) the second member of the variational equation is written as follows: where V } is defined as in (56); here | * | is the determinant of the Jacobian which toggles between the linear element  of curve ( 1 ) to real element .These formulas take into account the load spread over the entire surface and average load at the edges.

Validation Tests of Our Finite Element
The aim of this study is to investigate the accuracy of N-T thick shells theory for linear elastic shell by using Spherical Shell Finite Element (SSFE).We lay our investigation on a well-known benchmark as hemisphere under diametrically opposite charges given in Figure 1 used to evaluate the performance of a shell element.The computed deformed limit surface of a quarter of the hemisphere is shown in Figure 2 plotted with the Matlab R2015a tools and the displacement convergence results are shown in Figure 3 and Table 1.We monitor the displacement in load point  presented in Figure 1.
Here, the case study is that of a thin shell of hemisphere subjected to four opposite diametrically concentrated loads at the base proposed by Macneal and Harder [20] which  is discretized because of the symmetry of loads and the geometry.Both Kirchhoff-Love shell theory and N-T shell theory are computed using SSFE model and the respective results are analysed, compared with DKT12 and DKT18 proposed in [21] and SFE3 [22] then commented.Table 1 shows the displacement results at point  and Figure 3 perfectly describes their variation and the rate of convergence.The convergence properties of the method are clearly shown from Figure 3 and Table 1.Then SSFE converge as well as both the semifinite element (SFE) and Discrete Kirchhoff Triangle (DKT) elements for the membrane displacement at load point .

Scaling and Deviation. Scaling of the ratio 2𝜒 = ℎ/𝑅
given in Table 2 is proceeded on the range of the following values 0.006, 0.099, 0.12, 0.15, 0.175, and 0.2 of the spherical shell.Notice that the radius  is constant while the thickness varies with the ratio.We observe in Tables 1 and 2 and Figure 3 that, for the thickness ratio 0 ≺ 2 ≺ 0.099, the membrane displacement in load point  is the same for both K-L and N-T models.When the ratio 2 = ℎ/ is greater or equal to 0.099, the displacement computed for inextensible bending membrane from spherical equation of K-L and N-T is not the same.This means that, above 2 ≃ 1/10, both K-L and N-T approaches are different for all values of the thickness ratio.
We investigate now the deviation between K-L and N-T displacements in load point .With the variation of thickness ratio 2 = 0.06, 0.099, 0.12, 0.15, 0.175, and 0.2, the results plotted in Table 2 and Figures 4-7 clearly show that the deviation of displacement is encountered at the specified values of 2 above.This deviation increases with the number of meshes at the load point .We also observe that the deviation increases with thickness ratio.

Impact of the Transverse Stresses.
Surveys of various shear stresses in linear elastic thick shells can be found in the works of [3] where they are mathematically substantiated.We have obtained good convergence results for thin shells with the half side element  = 12 and a pressure in the load point   = 0.53 Pa for a spherical shell.For the thickness coordinate   ∈ [−ℎ/2, +ℎ/2], where ℎ = 1.2 m, the numerical results of the transverse shear stresses through the thickness are computed (see Table 3).After shear stresses investigation, we observe that they satisfy tangential stress-free boundary conditions at bottom to top surfaces of the panel.These results show predicting ability of SSFE finite element based on N-T's model.It also appears that the shear stresses vary through the thickness as compared to 3D thick shell equations.We can also appreciate that transverse shear stresses  13 and  23 are not negligible in thick shells.4.4.Discussion. 0 finite element discretization of the firstorder N-T model applied the PPR method to solve problems of gradient discontinuities across edges of triangular elements during interpolation.The displacement is two-degree polynomial of the thickness ratio  while the membrane strain tensor contains the change in the third fundamental form.
The convergence of SSFE has been clearly established and this element uses 12 degrees of freedom per triangular element which is robust and less greedy (in terms of memory) than DKT12, DKT18, and SFE that are based on the Kirchhoff-Love (thin shells) and Reissner-Mindlin (thick shells) approaches which neglect the third fundamental form in their shell kinematic equations [23].The constitutive law through the strain tensor contains the change of third fundamental form   as shown in the tables; the thickness ratio  impacts the behavior of the shell because of the significant energy contribution due to   as the ratio increases unlike the K-L classical model [3].
Recall that, in (17), we see that total deformation energy due to K-L and R-M models contains membrane deformation energy   and bending deformation energy   ; that is,  R-M =   +   .Equation (18) shows that the total deformation energy due to N-T model contains additional terms: coupled membrane and Gaussian bending (  ,   ) and Gaussian deformation energy   .Then  N-T =  R-M +   +   +   .
As we mentioned above, when the thickness ratio 2 is greater than 0.098, this additional energy (  ,   ) influences global deformation energy and shows the difference between N-T and R-M models applied to the spherical thick shells.
The investigation of the variation of the thickness ratio ℎ/ for certain values 0 ≺ 2 ≺ 0.099 proves that the Kirchhoff-Love, Reissner-Mindlin, and N-T classical models have the same contribution of total deformation energy.The additional terms containing the change of the third fundamental form do not have here any influence.The energies   and   disappear in the global deformation energy when 2 ≺ 0.1 or  2 ≺ 0.01 because they are inversely proportional, respectively, to 10 4 and 10 8 .In this case  N-T =  R-M .But if the thickness ratio is 0.099 ≤ ℎ/ ≺ 1, displacement results encountered for both models are different because   and   do not disappear in the global deformation energy.The consequence is that they enhance global deformation energy and improve the rigidity of the shell structure.Moreover, N-T's model brings us real facilitation to determine the distribution of transverse shear stresses through the thickness.

Conclusion
The finite element SSFE using the PPR method for the spherical shell described in this paper for  0 finite element triangularization has given accurate results and is skilled to design thicker shells.It has proved to be faster and uses less memory than other well-known methods used in the different benchmarks.The N-T model handles spherical thin and thick shell properly because it clearly shows how the change of the third fundamental form enhances the total deformation energy when the ratio  becomes greater.
Transverse stresses which have been predicted by the 2D governing equations of N-T model were calculated numerically and correctly.Structural engineers can therefore design more accurately spherical thick shells.

Figure 2 :
Figure 2: Deformed configuration of one quarter of the hemisphere.

Table 1 :
Table of results of displacement at point  for the hemisphere.

Table 2 :
Deviation of both Reissner-Mindlin (R-M) and N-T thick shell according to the scaling of the ratio ℎ/.

Table 3 :
Transverse shear stresses of a thick spherical shell under load point  convergence analysis.