An improved enhanced solid shell element for static and buckling analysis of shell structures

– This paper presents an improved higher order solid shell element for static and buckling analysis of laminated composite structures based on the enhanced assumed strain (EAS). The transverse shear strain is divided into two parts: the ﬁrst one is independent of the thickness coordinate and formulated by the assumed natural strain (ANS) method; the second part is an enhancing part, which ensures a quadratic distribution through the thickness. This allows removing the shear correction factors and improves the accuracy of transverse shear stresses. In addition, volumetric locking is completely avoided by using the optimal parameters in the EAS method. The formulated ﬁnite element is implemented to study the static and buckling behavior of shell structures and to investigate the inﬂuence of some parameters on the buckling load. Comparisons of numerical results with those extracted from literature show the acceptable performance of the developed element.


Introduction
In the last decades, the use of composite materials is becoming more widely used since they provide many advantages to structural designers.The stress and the strain fields in these structures are very complex; therefore, numerical methods have been developed to give an accurate prediction.In fact, in the case of plate and shell structures, to solve the shear locking problems in thin limit cases, the assumed natural strain method (ANS) have been used by several authors in the literature e.g.MacNeal [1], Bathe and Dvorkin [2] and Bucalem and Bathe [3] among many others.
Hauptmann et al. [4] assumed a linear distribution of the strain in thickness direction and developed a solid shell element formulation.The same assumption was taken by Sze and Yao [5] to formulate a locking-free solid shell element.
Another method was proposed by Simo and Rifai [6] named the enhanced assumed strain method (EAS).It allowed avoiding the thickness locking caused by coupling on the in-plane and transverse normal stress and normal strain responses [7].
Based on this method, Alves de Sousa et al. [8] proposed a volumetric and shear locking-free solid element a Corresponding author: mondherwali@yahoo.frformulation.The application of the EAS method for wide range of geometrically linear elastic structural analyses was presented in the works of Andelfinger and Ramm [9] and Rah et al. [10,11].The extension to non-linear aspects was carried out by Büchter et al. [12] and Bischoff and Ramm [13] among others.
The use of a single method cannot entirely solve the locking problems in solid shell element.Therefore, many studies have combined them for a better prediction of the structure behavior.Klinkel et al. [14] and Vu-Quoc and Tan [15] developed three-dimensional locking-free solid shell elements by using the EAS and ANS methods.Hajlaoui et al. [16] presented a solid shell finite element based on the classical ANS method to compute the shear locking effect and a nine parameters EAS method, which avoided completely the volumetric locking to study the bucking behavior of laminated composite plate with delamination.
In the literature, to the author's knowledge up to now, there is no solid shell formulation with high-order transverse shear enhancement.The only exception is the work of Quy and Matzenmiller [17], where the authors developed a finite element based on the theoretical foundations of high-order shear deformation theories.In this work, the authors decomposed the transverse shear strains into two parts: the first one is compatible and continuous, which can present a shear locking in thin limit structures, the second part is based on the EAS with parabolic function in terms of natural thickness coordinate.The membrane part in this work is not enhanced, which can present a numerical locking in the case of in-plane bending and in the nearly incompressible elasticity.
However, in the case of plates and classical shells, the higher order shear deformation theory (HOSDT) has been widely used by some authors e.g.[18][19][20] to improve the transverse shear strain and stress distributions.Recently, Wali et al. [21] present an efficiency three-dimensional double directors shell element for the analysis of functionally graded shell structures.The vanishing of transverse shear strains on top and bottom faces is considered in a discrete form inspired from the work of Dammak et al. [22].Thus, the third-order shear deformation plate theory (TSDT) is a particular case of the developed formulation.Nevertheless, the treatment of delamination with this HOSDT is complicated which is not the case with the solid shell elements.
This paper is an improvement of the solid shell formulation developed by Hajlaoui et al. [16].In the present formulation, the assumed natural transverse shear strain is enhanced with two parabolic functions in term of natural thickness coordinate as proposed in Quy and Matzenmiller [17].In addition, the enhancements of the membrane part and the transverse strain were used to avoid the in-plane bending and volumetric locking.
The remainder of this paper is organized as follows.Fundamental and finite element formulations are described in section two.In section three, numerical results and discussions of the finite element model are investigated in detail.Finally, some concluding remarks are presented in section four.

Solid shell finite element formulation
The developed solid shell element is an eight nodes hexahedral element with 3 degrees of freedom per node.Parameterizations of the solid shell material points are carried out in terms of curvilinear coordinates (ξ 1 , ξ 2 , ξ 3 ) = (ξ, η, ζ).The strain field is enhanced with the introduction of internal variables [6], as following: where E c and Ẽ are respectively the compatible part and the enhanced part of the Green-Lagrange strain tensor.

Compatible strains
The compatible part is arranged in (6 × 1) column matrix as follows: For the compatible part of the strain tensor, and to avoid shear locking, an assumed natural strain (ANS) method is used for the transverse shear strains E c 13 , E c 23 as proposed by Bathe and Dvorkin [2] where the transverse shear strains are evaluated at four mid-points of the element edges A = (-1,0,0), B = (0,-1,0), C = (1,0,0), D = (0,1,0), see Figure 1.
Also, an ANS method is used for the thickness strains E c 33 as used in references [14][15][16] where the transverse strain is evaluated using four collocation points defined in the reference surface : 1).
Then the compatible part of the Green-Lagrange strain tensors takes the following form: (5) where G ij = G i .G j and g ij = g i .gj are the metric coefficients of the reference and current configuration respectively.Here G i and g i are the covariant base vectors obtained by partial derivative of the position vectors with respect to convective coordinate (ξ 1 , ξ 2 , ξ 3 ) = (ξ, η, ζ) in reference and current configuration respectively.Matrix T , in Equation ( 5), is the transformation of the strain tensor from parametric coordinates to the local cartesian coordinates [14], written as See equation (6) next page.
where t ij = G i .T j and T j (j = 1, 2, 3) are a local orthonormal base vectors.The position vectors, within each element domain, in reference and current configurations are respectively where N is the tri-linear shape functions matrix, x n and X n are nodal coordinates.The displacement field, with the corresponding variation and increment, is interpolated in a same manner as follows where T is the vector of nodal displacements at the element level.Using Equation (5) and approximations (7) and ( 8), the virtual and Transverse shear strain points Transverse strain points incremental compatible Green Lagrange strain tensors are then given by where B is the strain interpolation matrix, relative to a node (I) and denoted B I is given by:

Enhanced strains
The enhanced part is related to the vector of the internal strain parameters α as: where Ẽ, δ Ẽ and Δ Ẽ are total, virtual and incremental enhanced Green Lagrange strain tensor respectively.The crucial assumption of the EAS method is the enforcement of the orthogonality conditions for the assumed stress field S and the enhanced strain Ẽ (Sect.2.3).This orthogonality conditions impose the following choice for the interpolation function matrix M to be expressed as follows M ξηζ dξdηdζ = 0 (12) where the subscript " 0 " means evaluation at the center of the element in the natural coordinates, is the Jacobian matrix.The interpolation matrix M ξηζ , in Equation ( 12), is expression in term of the parametric coordinates (ξ, η, ζ).The choice of matrix M ξηζ will be considered with 11 parameters Since the shear strain interpolation in thickness direction is parabolic, at least 2×2×3 Gauss quadrature integration rule is used.After including interpolation functions for enhanced strain fields, we obtain the solid shell element enhanced with 11 incompatible modes (C3D8C11).

Weak form and linearization
The variational framework of the enhanced assumed strain method (EAS), which is based on the three-field variational functional, in Lagrangean formulation, is given in references [14][15][16] among others.The expression of the variational functional, Π, is written as where ψ is the strain energy function and u, Ẽ and S are the independent tensorial quantities which are: displacement, enhanced assumed strain and assumed stress fields respectively.Vectors, F V and F S , in Equation ( 14) are the prescribed body force and surface traction respectively.When invoking the classical orthogonality condition, the number of independent variables in the original functional is reduced to just two u, Ẽ .The weak form of this modified reduced function may be obtained with the direction derivative leading to where S is the Piola-Kirchoff stress tensor given by: Equation ( 16) is a nonlinear equation that will be solved iteratively by the Newton-Raphson method which needs linearization.Using the enhanced approximation (11), linearization of Equation ( 16) is where L, H and K are given by with C = ∂ 2 ψ ∂E∂E the 6 × 6 three dimensional material tangent moduli, and K D is given by: In Equation (18), f int , f ext and h, are given by the following expressions K G is the geometric stiffness matrix Relative to a couple of nodes (I, J) The strain parameters Δα must be eliminated from Equation ( 18) at the element level, which leads to the element tangent operator, K T given by its classical partition into material and geometrical parts and the residual vector The basic equation of the buckling analysis is in the form of an eigenvalue problem where φ is the generalized global displacement eigenvector.This eigenvalue problem is solved using the subspace iteration method.

Numerical results and discussion
In this section we present several representative numerical simulations in order to illustrate efficiency of the proposed higher order shear strain enhanced solid shell element.All the numerical simulations have been performed with an in house nonlinear finite element program.
The performance of the proposed solid shell element is evaluated with several problems.The convergence of the results is compared to other well-known formulations.A listing of these elements, and the abbreviations used to identify them henceforth, are given in Table 1 (Appendix).

Analysis of stiffness matrix eigenvalue
In order to investigate the behaviour of solid shell elements in the incompressible range, an eigenvalue analysis for one cubic element with a side length 10 mm and Young modulus E = 1.0 MPa is performed.Since the material is assumed to be incompressible a Poisson coefficient ν = 0.499 is used.
Effective elements have a zero valued eigenvalue associated with each rigid body motion (the first six modes) and smaller ones for a particular deformation mode.To avoid the volumetric locking behaviour, it is important that the elements contain only one incompressible mode with a high valued eigenvalue (about 2500) for incompressible material.Table 2 shows the eigenvalues for 18 modes; the six zero eigenvalues for the six rigid body modes are omitted.
The elements Q8A3E5, Q8A3E7 and Q8S12 have undesirable high eigenvalues.The developed eleven parameters solid shell element C3D8C11 provides the correct eigenvalues.This proves the volumetric locking-free behaviour of the proposed C3D8C11 element that can be readily employed for the nearly incompressible analyses.

Cooks membrane problem
This linear example consists in a trapezoidal membrane clamped on one side while the other side is loaded (Fig. 2).
The panel has been analyzed assuming an elasticity modulus of E = 240.565MPa, a Poisson's ratio of ν = 0.499 and a thickness of h = 1 mm.The vertical displacements of the top edge node for various meshes and various element formulations are shown in Figure 2.
It can be observed from Figure 2, that the results obtained with C3D8C11 elements are considerably more accurate than those obtained with Q8A3E5, Q8A3E7 and Q8S12 elements.Our formulation exhibits much better accuracy for coarse meshes.The results obtained with Q8A3E5 and Q8A3E7 elements are nearly identical.Our developed formulation provides excellent results for the near-incompressible situations.This is due to the absence of volumetric locking.

Pinched hemispherical shell with 18 • hole
To investigate the capability of the developed elements to model the in-extensional bending and rigid-body motions and to overcome membrane and shear locking phenomena, we consider this classical pinched hemispherical shell.A hemispherical shell with an 18 • hole at the top is loaded at the free edge with two inward and two outward forces located 90 • apart (Fig. 3).
Material and geometric properties for this test are E = 6.825 × 10 7 MPa, ν = 0.49, radius R = 10 mm, and thickness h = 0.04 mm.Due to symmetry of the problem, only one quarter of the shell is modeled.The results obtained for the linear case are listed in Table 3 and shown in Figure 3.
The values for the displacements at the point of application of the force are compared to those obtained with Q8A3E5, Q8A3E7 and Q8S12 elements.It is well observed that the present results with C3D8C11 and those Top corner displacement (mm)  obtained with Q8A3E7 exhibits better accuracy even with coarse meshes.

Pinched cylinder with end diaphragms
To investigate the capability of the developed elements to model the in-extensional bending and complex membrane states of stress, we consider one of the more severe tests for shell and solid shell formulations.
A short cylinder, with two pinching vertical forces at the middle section, and two rigid diaphragms at the end, is modeled using one octant with appropriate symmetry boundary conditions (Fig. 4).The length of the cylinder is L = 600 mm, the radius is R = 300 mm, and the thickness is h = 3 mm.The material properties are : E = 3.0 × 10 6 MPa, ν = 0.3.The results are listed in Table 4 and shown in Figure 4.
In comparison to other solid shell formulations in the literature (Q8A3E5, Q8A3E7, and Q8S12), our developed formulation and those obtained with Vu-Quoc and Tan [15] show superior convergence behavior.However, the comparison of the present solid shell element with the one presented in Rah et al. [10] confirms that 7 EAS parameters are sufficient for wide geometrically linear elastic structural analyses.This proves the accuracy and efficiency of the present solid shell formulation for the linear elastic structural analysis of shells.

Buckling of isotropic plate
In this test, the performance of the method in incompressibility is examined.For this purpose an isotropic simply supported plate subjected to in-plane uniaxial compressive loads is considered (Fig. 5).
The material and geometric properties are: Young's modulus E = 1.0 × 10 7 MPa, edges a = b = 10 mm and thickness h = 0.1 mm.The results obtained for this case are listed in Table 5 and shown in Figure 5.
Figure 5 and Table 5 show the evolution of the normalized buckling load (N CriNum /N CriAnaly ) versus the Poisson's coefficient.In the case of thin plate, the critical analytic compressive force deduced from von-Karman plate theory is expressed as [23]: The results obtained with C3D8C11 elements are considerably more accurate than those obtained with Q8A3E5, Q8A3E7 and Q8S12 elements especially nearly the incompressible range (Fig. 5 and Table 5).This shows the performance of our element to avoid the volumetric locking     (Poisson-locking) in the near incompressibility limit and proves the accuracy and efficiency of the present solid shell formulation for the buckling analysis.

Buckling of laminate plate
The finite element model developed herein is validated by comparing the results with the three-dimensional elasticity solution given by Noor [24], Q8A3E5, Q8A3E7 and Q8S12 as shown in Figure 6, where a simply supported cross-ply [0/90/90/0] is considered, under uniform uniaxial loading.The material and geometric properties are: The present results with C3D8C11 have the closest results with 3D elasticity solution.They show more accurate values than those obtained with Q8A3E5 and Q8A3E7 elements.Comparing with results obtained with Q8S12 element our developed solid shell element C3D8C11 gives the same values.This proves the accuracy and efficiency of the present solid shell formulation for the buckling analysis of laminated structures.

Buckling of laminated plate with delamination
From the above results, we can conclude that the solid shell C3D8C11 developed in this paper gives the most accurate results comparing to other models.Thus, the C3D8C11 element will be used to investigate the buckling of laminated plate with delamination.
The following discussion extends the present analysis to study the effect of some parameters on the behavior of laminated composite plates with delaminations such as the stacking sequences, delamination size and aspect ratio.
To study the influence of the delamination size on the buckling behavior of a simply supported symmetric laminated plate, we consider a cross ply square laminated plate [0/90/90//0].It is subjected to a uniform in plan force with through-the-width delamination between the two last layers, as shown in Figure 7.The delamination in the through-the-width direction is performed with separated nodes.Moreover, due to the symmetry, only a quarter mesh configuration of the laminated plates with through-the-width delamination is studied.The material properties are: The buckling loads for different through-the-width delamination sizes of square plate with various values of (E 11 /E 22 ) are plotted in Figure 8.It is shown that the buckling load decreases as the delamination size increases.However, the buckling load tends to decrease slowly for small delamination (D/a ≤ 0.2), this phenomenon confirms that the presence of delamination reduces the strength of the laminate and the global buckling modes are converted into local buckling modes for small delamination size.
In this section, a four-layer cross-ply laminated plate subjected to in-plane uniaxial compressive loads, as  Figure 9 gives the buckling loads as a function of aspect ratio (a/b) and delamination size.It can be con-cluded that the bucking load decreases as the aspect ratio increases.
Figure 10 illustrates the variation of the buckling load versus the delamination size for symmetric cross ply laminated square plates.It is shown that the largest buckling load is obtained for the case of symmetric laminated [0/90/90//0], because the fiber orientation at the layer with the delamination is parallel to the in-plane loads.From the plots, it can be concluded that the influence of the delamination is more important for the symmetric laminate [0/90/90//0].

Conclusion
In this study an improved solid shell finite element based on the partition of transverse shear strain is developed.The first classical part, named ANS, is independent of the thickness coordinate and used to avoid the transverse shear locking in thin structures.The second is an enhancing part which ensures a quadratic distribution through the thickness.This later part allows removal of the shear correction factors and improves the accuracy of transverse shear stresses.The finite element formulation is completely free from the volumetric locking phenomenon, which occurs in the treatment of nearly incompressible elasticity.Numerical results show good agreement with analytical ones.The developed finite element model is validated by comparing numerical results with the threedimensional elasticity solution and other results from literature.The buckling behavior of composite laminates is investigated based on the use of solid shell element.Comparisons of numerical results with those from literature show the acceptable performance of the developed element.The influence of some parameters in the buckling behavior of a laminated composite structure with through-the-width delamination has been presented.The accuracy of the present results is due to a new optimal eleven-parameter solid shell element enhanced with the EAS method with high-order shear strain for the transverse shear strain linked to the ANS method.This formulation is being extended to study more complicated cases, especially those dealing with multiple delamination.

Fig. 2 .
Fig. 2. Cooks membrane problem; the vertical displacement of the top edge node.

Fig. 3 .
Fig. 3. Description and results of the pinched hemisphere with an 18 • hole.Symmetry is used and only one quadrant is modeled.

Fig. 4 .
Fig. 4. Description and results of the pinched cylinder with end diaphragms.Symmetry is used and only one eighth of the cylinder is modeled.

A.Fig. 5 .
Fig. 5. Description and results of critical buckling load of isotropic plate.

Table 1 .
Listing of solid shell elements.