Abstract

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 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 recovery method 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.

1. 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, and, as a consequence, the limit displacement reads , with , (where is the three-dimensional contravariant shell basis); the in-plane strain and the stresses depend on the in-plane strains and are computed from appropriate differential equations; see [3].

The form of the displacement clearly shows that 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) [46]; 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 -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) [810], 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) [1317], 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.

2. Materials and Methods

2.1. The N-T Model of Thick Shells

Let ( is midsurface, 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 and denote the covariant and contravariant basis of the midsurface and and , respectively, the covariant and contravariant basis of the shell. Thenwhere 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 while Latin indices take their values in the set . 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 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 and are the rotations angles. Strain tensors are miscalculated in some literatures because wrong basis vectors were used. The condition modifies the constitutive law which, expressed with the Lamé constants, now readswhere , or equivalently with Young’s modulus and Poisson’s coefficient , reads

Consider a shell of thickness , clamped on a part of its border , subject to volume forces and and to surface forces and on the rest of its border . Suppose the forces are sufficiently smooth; then the transverse stresses are solutions to the differential equations:where ; andwhere ; .

The generic point of the sphere is described bywhere and are curvilinear coordinates, , and ; 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 , be partitioned in two parts and the border of the shell with , and . We denote and . Suppose the shell is clamped on and subject to volume and surface forces as stated above; then the three-dimensional variational equation related to the equilibrium equation readswhere is the Sobolev space and: and denote, respectively, tensors and vectors scalar products. The variational formulation under Kirchhoff-Love’s approach which is given byis 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 . We should remember that the displacement calculated after isor in the Nzengwa-Tagne formatwhich suggests finite element implementation.

The linearized membrane strain tensor is expressed as follows:

Using the formulaswe can put (22) in the following form:wherewhere

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:

2.3. Resolution of the Transverse Stress Ordinary Differential Equation

Using expressions (27), (28), and (29), let , 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:

2.4. 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 . 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 finite element approximation of the solution , and letbe two triangular elements with a common node We denoteare the vector of nodal values in and , respectively. Let be the shape functions matrix; then restriction of in each element readswhere 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 , the set of polynomials of order ; then the gradient recovery PPR operator consists in defining such thatThis value also depends on the sampling points chosen. In the 2D case it is proved (cf. [14]) that if 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 , the existence and uniqueness solution are possible if ; 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 ; we denote ; we hereby determine so that .

Then, we obtainnotice that

For the transversal displacement, , the existence and uniqueness of a gradient recovery solution are possible if ; we fit quadratic polynômial using the same regular pattern:then, we scale by a factor with and , with respect to the six derivative values at the coordinates nodes of each element on the patch. Let , be the vectors of the matrix and consider the diagonal matrix and also the approximation vector of in a nodal point . Then we fit 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 , where . The union of the sampling points used in recovering the gradient at and the mesh nodes in the first layers around constitute the set of sampling point for recovering the gradient at ; see [15].

3. Discretization of N-T Model

3.1. Finite Element Space

The continuous truncated variational equation (18) is defined in the space . In order to define a finite element space, we begin by recalling the following results. Let be a triangle with nodes , , and . Let and be the sets of the first- and second-order polynomials with basis and , respectively, then and also generate and 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 bythe subspace of dimension (which is the number of nodes) and is the closure of , the subspace of dimension which is the number of unisolvent points (nodes and edge mid points).

The space is approximated with the finite element space of dimension

3.2. The Discrete Scheme

Let be a triangularization of the shell’s midsurface and be a triangle of vertices and midedges . Let as defined above. Then

By using the barycentric polynomials and writing , ; , , and ,

In the above formula, is the gradient operator and , , and are, respectively, the unknown values of in the midedges , , and . Let(i) and the coefficient of the approximation of the gradient when the vertex is internal of the mesh,(ii) and the coefficient of the approximation of the gradient when the vertex is in the boundary at the corner of the mesh,(iii) and 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 , we have

Using (43), (45), and (51) in which all three nodes of the triangular element are internal,  let  ; and we consider

Then the deformation vector can be written as where

3.3. Stiffness Matrix

Using local element, stiffness matrix can be expressed as follows:whereis the generalized behavior matrix.

Let and ; then we haveUsing (18) the second member of the variational equation is written as follows:where , and is defined as in (56); here is the determinant of the Jacobian which toggles between the linear element of curve to real element . These formulas take into account the load spread over the entire surface and average load at the edges.

4. 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 a standard test. This benchmark is usually used to verify the absence of membrane locking and good representation of rigid bodies motion. The hemisphere undergoes large rotations around the normal of the middle surface. Deformations of inextensible bending membrane are also important and this problem is therefore an excellent test to examine the ability of a shell element to represent the rigid and inextensible modes. The geometrical and mechanical characteristics are indicated for , the radius , the thickness , the angle subtended by the north pole of the hemisphere is , Young’s modulus is , Poisson’s ratio is , and the diametrically opposite loads at points and are . The limits conditions in the slib bound and the symmetric conditions are given by on the edge bound and on the edge bound .

4.1. Convergence

A reference solution presented in [20] provides for displacement in the direction of the loads as follows: . Only of the hemisphere 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 .

4.2. Scaling and Deviation

Scaling of the ratio 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 , the membrane displacement in load point is the same for both K-L and N-T models. When the ratio is greater or equal to , the displacement computed for inextensible bending membrane from spherical equation of K-L and N-T is not the same. This means that, above , 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 , 0.175, and , the results plotted in Table 2 and Figures 47 clearly show that the deviation of displacement is encountered at the specified values of above. This deviation increases with the number of meshes at the load point . We also observe that the deviation increases with thickness ratio.

4.3. 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 and a pressure in the load point for a spherical shell. For the thickness coordinate , 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 and are not negligible in thick shells.

4.4. Discussion

finite element discretization of the first-order 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, . 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 .

Where , , , and , which represents the portion of energy contribution, , and are constants which do not depend on .

As we mentioned above, when the thickness ratio is greater than , 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 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 or because they are inversely proportional, respectively, to and . In this case . But if the thickness ratio is , 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.

5. Conclusion

The finite element SSFE using the PPR method for the spherical shell described in this paper for 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.

Disclosure

Achille Germain Feumo is the corresponding author, Robert Nzengwa is co-first author, and Joseph Nkongho Anyi is co-second author.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding this research article.