Approximation of Linear Elastic Shells by Curved Triangular Finite Elements Based on Elastic Thick Shells Theory

We have developed a curved finite element for a cylindrical thick shell based on the thick shell equations established in 1999 by Nzengwa and Tagne (N-T). The displacement field of the shell is interpolated from nodal displacements only and strains assumption. Numerical results on a cylindrical thin shell are compared with those of other well-known benchmarks with satisfaction. Convergence is rapidly obtained with very few elements. A scaling was processed on the cylindrical thin shell by increasing the ratio χ = h/2R (half the thickness over the smallest radius in absolute value) and comparing results with those obtained with the classical Kirchhoff-Love thin shell theory; it appears that results diverge at 2χ = √1/10 = 0.316 because of the significant energy contribution of the change of the third fundamental form found in N-T model. This limit value of the thickness ratio which characterizes the limit between thin and thick cylindrical shells differs from the ratio 0.4 proposed by Leissa and 0.5 proposed by Narita and Leissa.


Introduction
In the field of structural mechanics the word shell refers to a spatial, curved structural member.The very high structural and architectural potential of shell structures is used in various fields of civil, architectural, mechanical, aeronautical, and marine engineering.The strength of the double-curved structure is efficiently and economically used, for example, to cover large areas without supporting columns [1].In addition to the mechanical advantages, the use of shell structures leads to aesthetic architectural appearance [2,3].Examples of shells used in civil and architectural engineering are shell roofs, liquid storage tanks, silos, cooling towers, containment shells of nuclear power plants, and arch dams.Piping systems, curved panels, pressure vessels, bottles, buckets, and parts of cars are examples of shells used in mechanical engineering [4].In aeronautical and marine engineering, shells are used in aircrafts, space crafts, missiles, ships, and submarines.Because of the spatial shape of the structure the behavior of shell structures is different from the behavior of beam and plate structures [4].
The considerable effort in the development of rigorous shell theories dates back to the early twentieth century [5][6][7].These shell theories reduce a basically three-dimensional problem to a two-dimensional one.Nevertheless, the analysis of shells with the aid of such theories involves complicated differential equations, which either cannot be solved at all [8][9][10], or whose solution requires the use of high-level mathematics unfamiliar to structural engineers.Therefore, many approximate shell theories have been developed, mainly on the assumption that the shell is thin, and to obtain generic analysis tools obviously some accuracy had to be traded for 2 Mathematical Problems in Engineering convenience and simplicity [4,11].Hence, it is not surprising that the development of the numerical formulations since the 1950s has led to a gradual cessation of attempts to find closed-form solutions to rigorous formulations [4].But, with today's availability of greatly increased computing power (also since the mid-twentieth century), completeness rather than simplicity is given more emphasis.
In computational mechanics, many numerical techniques are developed to solve governing equations of shells.We can identify for this purpose the generalized differential quadrature finite elements (DGQFEM) [11][12][13], discrete singular convolution method (DSC) [14], the finite difference (FD), the discrete volume finite element (DVFEM), and the finite elements method (FEM) under which solid-shells finite elements (SSFEM) and mixt finite element (MFEM) are found [15][16][17][18].Our approach is found under the MFEM category in which three-dimensional shell parameters are reduced to the midsurface of the shell of thickness ℎ.Some finite elements in this category describe well the behavior of shells under several loads.All these models of finite elements are link with either Kirchhoff-Love (K-L) or Reissner-Mindlin (R-M) hypothesis which neglected the effect of curvature in elastic stiffness [11,19].This is one of the main reasons why K-L and R-M theories of shells are incomplete and unsuitable to handle various types of shells [20].Structure engineers emphasis more on simplifying assumptions than limit analysis.Because of the third fundamental form that has been neglected in their theoretical approach of shells governing equations [21] and severe assumption on shell thickness, some mathematical terms without evident mechanic signification are introduced in the stiffness matrix in order to improve their efficiency [6,22].
After a brief presentation of Nzengwa and Tagne (N-T) kinematic equations of elastic thick shells, we develop a curved cylindrical 3-node finite element based on both K-L and N-T's shell models using strain tensor interpolation assumption [23].Next, through numerical implementation on classical benchmark, convergence of N-T's shell model is investigated.The deviation and the limit thickness ratio between the two theoretical approaches of thin shells and thick shells are established and discussed by applying successive scaling of a sample of thin shell.

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 geometrical or mechanical nature.The two-dimensional equations are deduced by applying asymptotic analysis on a family of variational equations obtained from an abstract scaled shell through multiple scaling of the initial three-dimensional equations.The theoretical elastic thick shells governing equations and limits displacements obtained from this approach (at the limit analysis) are more general because they contain additional terms to those found in Kirchhoff-Love model.The model also completes the thick shells theory of Reissner-Mindlin [24].
For a three-dimensional shell, the following bases ( 1 ,  2 ,  3 ), ( 1 ,  2 ,  3 ) are dual bases and ( 1 ,  2 ,  3 ), ( 1 ,  2 ,  3 ) defined on the midsurface also constitute dual bases of shell if  < 1.Let   ( 1 ,  2 , ) be the three-dimensional displacement field vector; it is defined as follows in the local coordinate system: The strain tensor for a three-dimensional shell reads From the limit analysis N-T demonstrated that the displacements field satisfies the equation  3 () = 0 and the unique solution is where (  ) are local membrane displacement components and  =  3 the transvers displacement: The general shell deformation tensor derived from the kinematics defined above is where   is the membrane strain tensor,   is the bending strain tensor, and   is the tensor of Gauss curvature.The strain tensor contains the change of the third fundamental form: In Kirchhoff-Love thin shells theory, the shell is assumed to be thin.Following that assumption, the three-dimensional displacements field vector   ( 1 ,  2 , ) is taken equal to the midsurface displacements field: It should be noted that the terms   = (        +   ∇  ) in the current model of thick shells disappear in the classical theory of thin shells because it is proportional to  2 and considered small.The thin shell theory, with the above assumption, cannot perform good shells behavior because the thickness has a considerable influence on the shells behavior.The thick shell model used in this survey appears to be suitable to handle both thin and thick shells because of its completeness.

Curved Shell Elements.
Curved shell elements are suitable to model the midsurface geometry more accurately.In the case of certain surfaces such as cylindrical shells, it means the exact description of the original surface.For more complicated cases, similar to the displacement field, the curvatures of the surface are approximated by interpolation functions.In this respect such elements belong to the parametric element types [25].

Triangular Curve Shell Element.
The triangular shell element is described in Figure 1; it can be double-curved (e.g., spherical structures) or single-curved element (cylindrical structures).In accordance with the N-T's equations of the theory of elastic thick shells, the geometrical properties of the cylindrical shell are the following: where  and  are curvilinear coordinates, 0 ≤  ≤  and 0 ≤  ≤ ,  ≡ cylinder radius, and , ,  are global coordinates.

Displacement Field.
Let   be the global displacement field vector defined as follows: With respect to the cylinder's geometry, the angles of rotations and the Gaussian curvatures are developed in function of membrane displacements  1 and  2 and the derivatives of the transverse displacement  , ,  , including curvature components    : which leads to where  1 ,  2 ,  1 , and  2 are geometrical properties of the midsurface of a shell within the respective coordinate-lines  and .
In the case of a cylindrical shell they are defined as follows [25]: The third component of rotation is added for a smooth transfer of element stiffness matrix from one system of coordinates to another.

Compatibility.
We recall the strain components in thick shell theory found in (3) as follows: The membrane deformation tensor   and bending tensor   are widely used in thin shell theory.But the tensor   is neglected [24].In classical thin shell and thick shell computation analysis, the contribution of this tensor (Gauss curvature tensor) to the deformation energy disappears in the stiffness matrix.As the shell becomes thicker, the contribution of Gauss curvature tensor in terms of energy is no more negligible as compared to that of the first two tensors used in classical thin shells.So this model is able to handle thick shells as well as the classical R-M models.
We express the strain components using the parameters of the cylindrical shell's midsurface.In detail, for the case of a cylindrical shell, we have 2.6.Interpolation.The approach we are using here is based on strain interpolation.This method consists in solving a system of differential equations from strain assumption.We have to find the displacements functions which perfectly capture the rigid-like motion; then particular displacements are calculated from approximation of the deformations components.These components are approximated such that membrane, shearing, and bending behavior are decoupled.So a pure bending state or a pure shearing state can be well represented.
In this section we consider  1 =  and  2 = V.
Since our displacement field is of six components, three rotations and three translations, we need six independent parameters to define our rigid body-like motion [18,25].
When ( 23) and ( 24) are taken into (18) and (19), it comes that the rigid body-like motion is well handled by the following displacement components: The displacements found in (25) are expressed in matrix form: where 0 is the vector of unknown coefficient which captures the rigid body-like motion: 0 is the matrix of interpolation functions.
The general solution displacements field  is the sum of the displacements field of rigid body-like motion and that of a particular solution of the deformation   : where The strain components are interpolated in such a way that interference between shear strains and bending is not possible.The following displacement components satisfy strains defined above:

= (
0 0 0 0 0 0 0 0 0 0 0 0   0 0 0 0 0 0 0 ) . ( is the interpolation functions matrix related to the particular solution of displacements field.From the total displacement vector field expression above, we have 18 unknowns.It then requires 18 displacement parameters to find all the above unknowns (see Figure 2).So the total displacement vector over the triangular element is In the local (curvilinear) system of coordinates, each node has six (06) components of displacement.

Stiffness Matrix.
In order to calculate the stiffness matrix, we have to formulate the variational problem over a domain.
Let  be the border of the domain; then let  =  0 ∪  1 be partitioned in two and the border of the shell Ω = Γ 0 ∪ Γ 1 with Γ 0 =  0 × {−ℎ/2, ℎ/2}, and We define Γ − =  × {−ℎ/2} and Γ + =  × {ℎ/2}.We suppose here that the shell is clamped on Γ 0 and loaded by volume and surface forces as stated above; the three-dimensional variational equation related to the equilibrium is is volume forces in the domain Ω and  is surface forces in the domain Ω.
The constitutive law of the linear elastic homogenous material is where  =   =     + (    +     ),  and  are Lamé coefficients which depend on intrinsic properties of materials.
Replacing (34) in (33), the problem is stated as follows: The strain tensor is defined at (7), if substitution is done in the left side of the above equation, it can now be written: Naturally  3 () = 0 at the limit analysis.Then at the midsurface,   = ( −1 )   ( −1 )     .By using Taylor's expansion on and truncating at  = 1, we obtain the best first-order twodimensional Remark 1.For Taylor's expansion used above, if we truncate (37) for  = 0, the bilinear continuous form (, V) will take the value: which is found in the literature to be K-L shell deformation energy.It is included in  1 (, V) which is the N-T shell deformation energy.In what follows, we compute both  N-T  and  R-M  which are, respectively, stiffness matrix related to  1 (, V) and  0 (, V).
For the computing aim,  1 (, V) is transformed For this purpose, components of strain tensor are expressed as a product of the matrix of interpolation functions' gradient and that of the total displacement vector over the domain.From (39), the stress tensor components are deduced in function of the product above.
Also, the bending stress and strain are expressed: Finally as previous, the Gauss curvature tensor is expressed as Total displacement is stated as the sum of particular solution and homogeneous solution of our system of strain differential equations.
where  = [ 0   ] is the matrix of interpolation functions and The left-hand side of the variational problem reads This is written as follows: where  (⋅)  is the global stiffness matrix.The midsurface is symmetric in shell thickness: that is, the midsurface is at +ℎ/2 from the top surface and at −ℎ/2 from the bottom one.(49) The right-hand side of the best first-order two-dimensional variational equation reads where  and  are defined as in [24] and  is the force vector from distributed load   + .
By equating  1 (, V) to (V), we obtain the structural equilibrium equation over the whole area as follows: Over a single-curved triangular element of area   the elementary stiffness matrix   and elementary force vectors   can be locally implemented.
By regular assembling process of elementary stiffness matrix   and elementary force vectors   , we then expressed the global stiffness matrix and force vector over the whole area as follows: is total number of elements.

Numerical Study
A numerical study which aims to quantitatively investigate the accuracy of different thick shell theories for the linear elastic analysis of cylindrical structures and to provide a basis for the long-term analysis is presented.
One of the problems frequently used to evaluate the performances of a shell element is that of the cylindrical roof (of Scordelis-lo) [4,26] subjected to its self-weight described in Figure 3.The flat rims are free and the curved edges rest on rigid diaphragms in their plans.The geometrical and mechanical characteristics are indicated for ℎ/ = 0.01, /ℎ = 200.The length of the cylinder  = 600 cm; its radius is  = 300 cm; the angle subtended by the roof is 2 = 80 deg; the thickness is ℎ = 3 cm; Young's modulus is  = 3.10 + 10 Pa and Poisson's ratio is ] = 0.

3.1.
Convergence.The quarter of the roof is discretized by considering triangular single-curved elements with  = 4, 6, . . .elements on edges  and .The transverse displacements at  and  are plotted (Figure 4) and compared with other models of finite elements: see Figure 5.
The computed deformed limit surface is shown in Figure 4. Here, displacement-convergence results are shown in Figure 5 and Table 1.In the diaphragm case, we monitor the displacements under the self-weight at points  and .The convergence properties of the method are evident from Figure 5 and Table 1.As may be seen, the CSFE3 converges as well as both the semifinite elements (SFE) and the DKT elements for the displacements under the loads.

Scaling and Variation of Displacements.
A successive scaling is carried out on the thickness ratio ℎ/ over the Scordelis-lo roof.We have used the following range of ratios: 0.1; 0.2; 0.3; 0.325; 0.4; 0.5, because it certainly belongs to both thin shells and thick shells ranges of thickness.The radius  is constant while the thickness ℎ varies with the ratio.The scaled shell obtained here has been implemented using both K-L theory and that proposed by N-T.The results are plotted for each ratio as shown in Figure 6.
We observe from the diagrams shown in Figure 6 that for the thickness ratios ℎ/ = 0.1, 0.2, and 0.3, the transverse displacements at point  are the same for single-curved triangular computation of both K-L and N-T models (see Figures 6(a), 6(b), and 6(c)).One can verify that, at any other points of the roof structure, these same observations are satisfied.For thickness ratios ℎ/ greater or equal to 0.325, we observe that transverse displacements computed from the shells equations of K-L and N-T are not the same (see Figures 6(d), 6(e), and 6(f)).We can still investigate the exact value ℎ/ in ]3/10 13/40[ from which displacements results issued from these two approaches of shells theories diverge in a cylindrical wall and probably other shells.

Deviations
. N-T's equations for elastic thick shell are more general [24] and with regard to that, we assume the displacements issued from this model of thick shell to be more accurate.We can then investigate the deviations encountered when using K-L model to evaluate a shell.The investigation consists in calculating the deviations (difference between K-L displacement and N-T displacement) for each mesh-step (number of elements per half side) with respect to the series of thickness ratios ℎ/ = 0.1; 0.2; 0.3; 0.325; 0.4; and 0.5 at the same point.This investigation was carried out on two different points  and  (see Figure 3) and the results are plotted in Figures 7(a From diagrams shown in Figures 7(a) and 7(b), we observe that (i) at points  and  of our cylindrical roof, deviations are encountered from thickness ratios ℎ/ above 0.3; (ii) deviations increase with mesh-steps at point  and decrease at point .

Discussion
We have investigated in this framework the influence of the third fundamental form proportional to  2 (with  2 = (ℎ/2) 2 ) included in N-T's model of elastic thick shells.This model was implemented using cylindrical shell triangular finite elements (CSFE3).Based on the results obtained, the following comments can be made.
The convergence studies of CSFE3 show that this curved finite element of three nodes only converges as well as SFE, DKT12, and DKT18 which are robust and greedy (memory) [27].The last three finite elements models are based on theoretical approach of Kirchhoff-Love (thin shells) or Reissner-Mindlin (Thick shells) which neglect the third fundamental form in their shell kinematic equations.In addition to that, these classical shell theories are based on two powerful assumptions of K-L or R-M [28,29].
The theoretical equations of thick and thin shells issued from these assumptions can be questioned because the transverse fiber (normal to the reference surface) has been forced to follow a specific behavior without any scientific or technical background that justifies their kinematic equations [20,24,26].Consequently some quantities without evident mechanical significations are added in the formulation of finite elements models (like DKT, SFE . ..) in order to improve computational accuracy [7,26].
In the contrary, our finite element (CSFE3) is based on N-T's theoretical model of thick shells which is deduced from  3 () = 0 obtained at limits analysis of 3-dimensional shell equations [24].Resultant kinematic equations of shell here are more general since they contain additional terms to those found in classical equations of shells.They are mathematically and technically justified without ad hoc assumption on   transverse fiber's behavior.This can explain why a simple 3-node curved triangular finite element obtained efficient convergent result.The investigation with a 6-node curved triangle will surely improve the convergence rate [7,26].The successive scaling over the benchmark structure has presented divergence of displacements results as ℎ/ increases for K-L kinematic model of shells and that of N-T.This is due to the influence of the additional terms found in N-T's shell theory in the deformation energy of the structure.For homogeneous and isotropic elastic shells with a midsurface between −ℎ/2 and +ℎ/2, the total deformation energy contains membrane stiffness   , bending stiffness   , coupled stiffnesses   and   (membrane and Gaussian bending), and finally Gaussian stiffness   .
For Reissner-Mindlin and K-L classical model of shells, the global stiffness matrix  R-M  contains only   and   while in N-T's model, the global stiffness matrix  N-T  contains additional terms (  +   +   ) to those found in  R-M  .The proportion of energy contributed by each stiffness matrix to the global deformation energy is as follows:   =   (ℎ/);   =   (ℎ/)(ℎ/) 2 ;   =   (ℎ/)10 −1 (ℎ/) 2 ; and   =   (ℎ/)10 −1 (ℎ/) 4 .We monitored that for thickness ratio ℎ/ less or equal to 0.3, that is, (2) 2 < 1/10, the couple stiffness   and Gaussian stiffness   disappear in the global deformation energy because they are inversely proportional, respectively, to 1000 and 10.000.We then obtain  N-T  ≈  R-M  .This explains why deviations are zero for thickness ratio less than 0.3 (see Figure 7).For thickness ratio greater or equal to √1/10 displacements results encountered for the two models are different because the coupled stiffness   and the Gaussian stiffness   which are now inversely proportional, respectively, to 100 and 1000 cannot disappear anymore in the global deformation energy of the structure.This additive proportion of energy into the global deformation energy found in this model improves total stiffness and increases the rigidity of the structure.

Conclusion
A finite element model for predicting the elastic thick shell behavior of cylindrical structure has been developed in this work.It has been quantitatively shown that the thickness ratio (thickness on radius) plays important role in their elastic behavior and structural safety.Hence, a reliable control of this parameter is required.N-T's theoretical approach for the modelling of the displacement and stains in cylindrical structures have been examined and compared with that of K-L for the case of self-weight loading cylindrical roof.
The comparison has revealed that N-T's shell theory is suitable for the analysis of cylindrical thin and thick structures from the perspectives of both accuracy and simplicity.A 3-node cylindrical shell finite element has been developed for the simulation of static behaviour of shells in general, which requires the ability of well handling thickness ratios from the vicinity of zero to one.
The ability of the CSFE3 model to deal with a wide range of ratios has been demonstrated through numerical examples.
It has been found that the use of K-L theoretical approach based on those assumptions leads to inaccurate results as the thickness becomes greater than the square root of 0.1.
It has also been found that above the value square root 0.1, the deviations linearly increase with the increase of the thickness ratio.
The issues addressed in this work highlight the need for N-T's approach to be implemented for the elastic shells analysis of cylindrical structures.
The analytical model and the numerical results presented here contribute to the establishment of a foundation of theoretical knowledge required for reliable analysis, effective design, and safe use of cylindrical walled structures.
Finally, the paper provides a platform for double-curved finite elements model to be developed for studying the influence of thickness ratio on all shell geometries.

Figure 3 :
Figure 3: Definition of Scordelis-Lo problem.Cylindrical roof subjected to its self-weight.

Figure 4 :
Figure 4: The undeformed and deformed configurations of onequarter of the roof.

Figure 7 :
Figure 7: (a) Evolution of displacements deviations at point .(b) Evolution of displacements deviations at point .

Table 1 :
Comparison of transverse displacements at  and .