An 8-Node Shell Element for Nonlinear Analysis of Shells Using the Refined Combination of Membrane and Shear Interpolation Functions

An improved 8-node shell finite element applicable for the geometrically linear and nonlinear analyses of plates and shells is presented. Based on previous first-order shear deformation theory, the finite element model is further improved by the combined use of assumed natural strains and different sets of collocation points for the interpolation of the different strain components. The influence of the shell element with various conditions such as locations, number of enhanced membranes, and shear interpolation is also identified. By using assumed natural strain method with proper interpolation functions, the present shell element generates neither membrane nor shear locking behavior even when full integration is used in the formulation. Furthermore, to characterize the efficiency of these modifications of the 8-node shell finite elements, numerical studies are carried out for the geometrically linear and non-linear analysis of plates and shells. In comparison to some other shell elements, numerical examples for the methodology indicate that the modified element described locking-free behavior and better performance. More specifically, the numerical examples of annular plate presented herein show good validity, efficiency, and accuracy to the developed nonlinear shell element.


Introduction
The 8-node isoparametric serendipity shell finite elements were suffering from locking effects due to smaller thickness.To avoid this deficiency, reduced and selective integration techniques in the finite element method have been proposed; however, spurious zero-energy kinematic modes occur and may disturb the finite element response in a mesh.
Hinton and Huang [1] developed the 8-node element model, but the element showed poor accuracy rather than the 9-node quadrilateral element.Lakshminarayana and Kailash [2] presented an 8-node shell element with locking free.In order to describe the locking problem, they used appropriately chosen interpolation functions based on Hinton and Huang's concept.
Bathe and Dvorkin [3] proposed an 8-node shell element-MITC8 to avoid membrane and shear locking problem.The strain tensor was expressed in terms of the covariant components and contravariant base vectors.The performance of this element was quite satisfying and suggested as the promising results in very complex shell structures.Bucalem and Bathe [4] had improved in previous publications the MITC8 shell elements [3], and the results provided conservative and unconservative performance in the finite element method.
Kim and Park [5] and Kim et al. [6] presented an 8-node shell finite element.In an 8-node shell element, the persistence of locking problems was found to continue through numerical experiments on the standard test problem of MacNeal and Harder [7].In recent year, Han et al. [8] studied the new combination of sampling points for the assumed natural strain and concluded that while it performed quite effectively in some cases, in a few analyses of very thinwalled structures, the accuracy of the element showed less than the predicted results.

Mathematical Problems in Engineering
Han et al. [8] used the membrane strain sampling points by Bathe et al. [9] and Han et al. [10].The shear strain sampling points by Hwang [11] are selected in Han et al. [8].
In present study, the interpolation functions by Polit et al. [12] for the sampling points (Han et al., [8]) are used and the strain component of center point is replaced by mean of the components at two points (Bathe and Dvorkin, [3]).
The primary goal of this paper is to propose an improvement of the most useful curved quadrilateral shell finite element, which is clearly, from a practical point of view, the 8-node element and then in order to improve the 8-node ANS shell element, a new combination of sampling points and shaper functions are adopted for the analysis of plates and shells.Also, in order of validation of the present shell element models, the numerical examples are investigated and compared with those solutions from the previous literatures.

Improvement of Shell Element
2.1.Kinematics of Shell. Figure 1 shows the geometry and current kinematics of an 8-noded shell element with six degrees of freedom.
Kinematic equations for the first-order shear deformation theory including extension of the normal line can be obtained from the 3D equations of the theory of elasticity by using a well-known first-order approximation of a vector function with respect to the coordinate  3 (Rikards et al, [14]).Further well-known geometric relations for the shell in normal coordinates are also used.
Let us assume that vector P characterizes the position of an arbitrary point of the shell in the initial reference state (see point  in Figure 2) and vector Q is the position of the same point (point   ) in the deformed state.The position of a point at the mid-surface of the shell (point ) in the initial state is characterized by vector P, and position of the same point in the deformed state (point   ) is characterized by vector Q.Normal curvilinear coordinates   = [  ,  3 ] at the midsurface of the shell in the initial state are defined by the righthanded triad of the base vectors [a  , a 3 ].Here, a unit vector a 3 is normal to the mid-surface of shell.Therefore (see Figure 2), Similarly, curvilinear coordinates deformed state are defined by the triad of vectors [A  , A 3 ].In the deformed state, the vector A 3 may not be perpendicular to the mid-surface of shell.
Vector function Q can be expanded in a Taylor's series with respect to coordinate  3 normal to the mid-surface of shell as follows: Here, ∇ ⊗ Q = A  ⊗ a  = G is a second-order tensor (see Luri' e, [18]), which characterizes the gradient of strains with respect to coordinate  3 , where ∇ is a Hamiltonian operator and ⊗ is a dyadic product of the tensors.Further, the notations for displacement at the mid-surface u and displacement of an arbitrary point of shell u are introduced.The following expressions can be written (see Figure 2) as From expressions (1)-( 3), the representation of the displacement u of an arbitrary point of the shell for the first-order approximation can be obtained from where  is vector of rotation at the mid-surface of shell where The incremental form of the displacement field may be written in terms of the nodal incremental vector Δu  as where } and H  is the rotational matrix of the normal vector at node where T  R is the transformation matrix between the initial shell normal and the deformed normal can be written as the result of a sequence of finite rotations  1 ,  2 , and  3 as follows: The incremental membrane and bending and transverse shear strains can be separated into linear and nonlinear parts (see Han et al., [10]).The linear strain parts can be expressed as: The shell theory outlined above is the so-called first-order shear deformation theory with six degrees of freedom.

Strain Energy and Stiffness of Shell.
The strain energy  of the shell represented as a three-dimensional body is given by the expression, where in curvilinear coordinates, the contravariant stress tensor   is contracted with the covariant strain tensor where  is a volume element of the shell.The Hooke's law for each layer can be written as Here,   is the component of the elasticity tensor in the mid-surface metrics.Substituting ( 12) into (11) as the strain energy  can be expressed as The plane stress condition is enforced.For single layer plates and shells, these equations are quite straightforward.The membrane forces, the bending moments, and the transverse shear forces can be obtained by integrating the relevant stresses through the thickness using the equivalent constitutive equations: where The membrane-bending D  and shear strain rigidity D  have 6 × 6 and 2 × 2 matrices, respectively where C 1 is the elastic constitutive coefficient where   is the shear correction factor (  = 5/6) and  is the shear modulus.When the entire cross section is isotropic elastic, then the rigidity matrix reduces to the following form: where C 1 and C 2 are the elastic constitutive coefficients of membrane-bending and transverse shear parts.
In this study, for the improved efficient 8 node shell element, the usual 8-nodes of Lagrangian displacement interpolations are employed, and the various combinations of assumed natural strain interpolation functions are used.Figure 3 lists various patterns of sampling points that can be employed for membrane and in-plane shear and transverse shear strain interpolations for the improved 8-node shell element.Based on Figure 3, the  pattern is used for membrane ( and ), and the  pattern is used for membrane (), as well as transverse shear ().The  pattern and  pattern are used for in-plane and transverse shear, respectively.Han et al. [8] used the transverse shear interpolation functions by Huang [11] in the  patterns.The interpolation functions by Polit et al. [12] are used in the  * and  * 6 patterns.In the  6 and  * 6 patterns, the strain component of center point is replaced by the mean of the components at points  1 and  2 (Bathe and Dvorkin [3]).The six cases of the combinations of various sampling points are used in the analysis.

Incremental Equation of Equilibrium
The generalized Hook's law at large strain does not represent an approximate material behavior description because stressstrain relation is non-linear.From the practical point of view, Hook's law is only applicable to small strain, for which constitutive tensor is constant coefficient.
The following incremental equilibrium equation is obtained where superscript  which is generally used as the current configuration is ignored in (18) and superscript  + Δ is the adjust incremented configuration; +Δ   ext is the external virtual work in  + Δ.
The total tangent stiffness comprises the material stiffness and the geometric stiffness.The linear part of the Green strain tensor is used to derive the material stiffness matrix, and the non-linear part of the Green strain tensor is used to derive the geometric stiffness matrix.

Linear Element Stiffness Matrix.
If the strain displacement ( 8) is substituted into (18), the linearized element material stiffness matrix (K L ) is obtained as follows: Mathematical Problems in Engineering 5 Finally, the element stiffness matrix has 6 × 6 size on the reference surface of shell element where the submatrix of [K L ] is shown in Han et al. [10].

Geometric Stiffness Matrix.
In order to obtain an accurate geometric stiffness matrix, the stresses should be evaluated accurately.The accuracy of the computation of stresses for formulation of geometric stiffness matrix is maintained by obtaining the same interpolated strains in the computation of linear stiffness matrix.The stresses are computed at the integration points based on these strains.Substituting the non-linear part of strain into (18), the following geometric stiffness matrix ([K G ]) is obtained: The geometric stiffness matrix in the natural coordinate is analytically integrated through the thickness.By the transformation of the natural to the global frame, the element geometric stiffness matrix is obtained on the global frame with 6 × 6 submatrix as follows: where, the sub-matrix of [K G ] is shown in Han et al. [10].
Then the final assembled incremental non-linear equilibrium equation can be written as where, F and F are the external and internal forces, respectively.
The equilibrium equation must be satisfied throughout the complete history of loading, and the non-linear processing will be stopped only when the out of balance forces are negligible within a certain convergence limit.If it is necessary to extend the stability analysis beyond the limit point, that is, in the so-called post-buckling range, appropriate solution procedures must be applied.One approach is to use the arc-length control method in conjunction with the Newton-Raphson method to extend the stability analysis beyond the limit point, by Crisfield [21].

Torsional Effect.
The element stiffness matrix may be written in a matrix form using the equivalent constitutive equations.Finally, the element stiffness matrix has 6×6 size on the reference surface of shell element.The torsional stiffness term was formed as described by Kanok-Nukulchai [22] and added to the stiffness term.
In this study, based on the procedure proposed by Kanok-Nukulchai [22], the drilling degree of freedom will be tied to the in-plane twist by a penalty functional through additional artificial strain energy as where   is a parameter to be determined (the value of 0.1 suggested);  is the shear modulus;   is the volume of the element; and  is the volume element.After integration throughout the thickness, ( 24) can be written as If   ℎ is chosen to be large relative to the factor ℎ 3 (which appears in the bending energy), (25) will play the role of penalty function and result in the desired constraint: At the Gauss points, A two-by-two Gauss integration scheme is applied for the evaluation of the torsional stiffness in order to avoid the overconstrained situation.To derive a torsional stiffness from (25), the local variables are expressed in terms of global nodal variables, u, by shape functions.This gives (25) in the form (Kanok-Nukulchai [22]) The torsional stiffness term K tL is added as described by Kanok-Nukulchai [22].

Numerical Results of the 8-Node Shell Element
In this section, we describe numerical tests of improved 8node shell elements listed in

Patch Test.
In the study, the basic patch tests proposed by Simo et al. [13] were performed, and the results illustrated in Figure 4.A patch of five elements is used in this study.Also, Boundary and loading conditions are illustrated in Figures 5 and 6, respectively.Tables 2, 3, and 4 present the normalized solutions of nodal displacements on the right edges.  = / = 0.12 × 10 −4 ).

Combinations of various sampling points
*  6  * 6 Normalized solutions 1.000 1.000 1.000 1.000 1.000 1.000The normalized solutions are presented in the nondimensional form using the following equation: Patch test results indicate that various sampling points for assumed natural strain method in Tables 2 to 4 can represent fields of constant moment, constant in-plane tension, and transverse shear forces.In some cases, the present shell element can pass the patch test when patterns , ,  6 , and  * 6 are used in all cases.In other cases, the shell element cannot pass the patch test when pattern  is used in pure transverse shear case.Therefore, we found that it is very important to apply the best combination of sampling points for 8-node shell elements.

Bending of Rectangular Plate.
The simply supported and clamped rectangular plate problem under uniform and central point loadings is applied to test shear locking by changing aspect ratios.Two aspect ratios of / = 1 and 5 were considered, and a quarter was modeled due to symmetry.To test the general applicability, the plate is analyzed using both a rectangular mesh and a distorted mesh.Reference solutions (MacNeal and Harder [7]) are given in Tables 5-10.The reference vertical deflection at the center of the simply supported plate (shown in Figure 7) under uniform load is 4.062, and the clamped plate under concentrated load is 5.60.In Tables 5 and 6, the numerical results obtained by using different types of existing elements are listed as follows.
(b) For the Case of / = 5.0.
The reference deflection at the center of the simply supported plate under uniform load is 12.97, and the clamped  plate under concentrated load is 7.23.In Tables 7 and 8, the numerical results obtained by using different types of elements are listed as follows.
Locking problems occurred when the  pattern is applied to transverse shear strain in thin plates (Han et al. [8]).This particular example showed the limit of 8-node shell element   square plate, /ℎ = 20000, is analyzed using both a rectangular mesh and a distorted mesh.A clamped boundary condition is chosen because it is considered to be more severe compared with simple supports.Two types of loading are considered, concentrated load and distributed load.In Figure 8, the problem geometry and material properties are provided.A 4 × 4 mesh is used on one-quarter of the plate.The results are compared in Table 9 to those provided by Timoshenko and Woinowosky-Krieger [25].The accuracy of the results obtained for both cases is maintained and compared with Choi et al. [20] and Kim et al. [16].
The reference vertical deflection at the center of the plate is the simply supported plate under uniform load which is 4.062 and the clamped plate under concentrated load which is 5.60.In Tables 9 and 10, the numerical results obtained by using different types of existing elements are listed as follows.
The cases of  and  * 6 are compared for the test of distorted mesh.It is also found that the locking phenomenon does not happen and the solutions with rectangular mesh are of similar accuracy.

Curved Beam Problem.
The cantilever beam is clamped at one end and loaded by a unit force at the other.The force is applied either along the in-plane axis or along the out-ofplane axis as seen in Figure 9.In the curved cantilever beam, combination of the principal deformation modes is evoked by a single in-plane or out-of-plane shear load at the free end.Note also that the element shape is quite rectangular, which will test the effect of slight irregularity.Theoretical displacement of free end under in-plane load is proposed 0.08734 by MacNeal and Harder [7].Recently, a different solution of 0.08854 was given by Young.[26].Reference solutions (Young [26] and MacNeal and Harder [7]) are given in .
Results presented in Tables 11 and 12 show the outstanding performance of the proposed element.

Twisted Beam Problem.
The initial geometry of the beam is twisted, but the initial strain is equal to zero.The twisted beam problem (shown in Figure 10), which was introduced by MacNeal and Harder [7], was proposed to test the effect of element warping.The performance of elements of thickness 0.32 was investigated under in-plane and out-of-plane shear loads.The warp of each element is 15 ∘ .Numerical results in Tables 13 and 14 are listed with displacements of free end normalized to the reference solution (MacNeal and Harder [7]).
Results are presented in Tables 13 and 14, showing more correct results than the result of references in spite of the 6×1 mesh.

Pinched Hemispherical Shell with 18 ∘
Hole.The hemispherical shell with an 18 ∘ hole is loaded by two pairs of equal but opposite external forces; see Figure 11.The shell undergoes an almost inextensional deformation, and a membrane locking can destroy the solution.For the shown FE mesh, the elements are flat and trapezoidal.Only a quarter of the shell is modeled because of the double symmetry.
For this problem, the reference deflection at point A in the -direction is 0.093 (Simo et al. [13]).The geometry and material properties are shown in Figure 11.The numerical results are given in Table 15.
This indicates that the results obtained by our elements are more accurate than the solution by the QUAD8 * * element.Note that the error of 8 × 8 mesh for the QUAD8 element is greater than 10%.12, is yet another challenging test for the ability of the element in modeling inextensional bending modes and complex membrane states.Most of the 4-node shell elements do not show good behavior in this test, but the present 8-node solid element and others passed this test.The vertical deflection of point A was investigated.The reference value is 0.18248 × 10 −4 (Belytschko et al. [27]).The results with different meshes and elements are presented in Table 16.A good performance of the current elements was again observed.The cases of  and  * show that the results are more conservative than the reference solution.However, the solution in the case of  * 6 is more accurate than that in the reference element XSHELL-8-ANS.

Nonlinear Analysis of Slit Annular Plate.
To validate the present shell element for geometrically nonlinear analysis, a numerical example of annular plate under end shear force is solved.We examine an annular clamped plate subjected to a distributed transverse shear force (see Figure 13).This good benchmark problem was considered by Buechter and Ramm [28], Brank et al. [29], and Sansour and Kollmann [30], among others, and Sze et al. [31], presented ABAQUS results for the problem.In particular, we use the ABAQUS results presented in a recent paper of Sze et al. [31] for popular benchmark problems of nonlinear shell analysis (because of the tabulated data).In Figure 13, the geometry and elastic material properties for the isotropic case are given.
The line force  is applied at one end of the slit, while the other end of the slit is fully clamped.In the present example, the plate is modeled by polar coordinates.The problem, however, can be solved by using the standard Cartesian coordinates as well.A regular mesh of 5 × 40 elements is considered in the present analysis.Solutions obtained by the mesh of 2 × 16 are too stiff and suffer from locking.The arclength method is used to obtain the load-deflection curve in  these examples.The shear load versus displacement curves for three characteristic directions are depicted in Figure 14.Solutions obtained with the present formulation are in good agreement with those obtained by Sze et al. [31].Figure 15 shows the deformed configuration of the isotropic annular plate for  = 1.45.

Nonlinear Analysis of Hinged
Shell.This is a widely used example for testing geometrically nonlinear shell elements.In this study, a quadrant of the shell is modeled with 12 × 12 mesh sizes, as shown in Figure 16.Three different thickness values are considered; that is, ℎ = 3.175, 6.35, and 12.7 mm, respectively.The problem has been considered in Sabir and Lock [32], Horrigmoe and Bergan [33], Rhiu and Lee [34], Sze et al. [31], and Arciniega and Reddy [35].The arc-length control method is used to trace the equilibrium path, and a tolerance coefficient of 1.0 × 10 −6 is used.The total length is  = 508 mm, curvature of the circular edges (1/) is 1/2,540 mm, and the circumferential length is 508 mm.The material properties are Young's modulus  = 3, 102.75 N/mm 2 , and Poisson's ratio V = 0.3.
For the case of ℎ = 12.7 mm, the curve of the central deflection versus load is given in Figure 17, showing the good correlation between the present solution and the existing solutions.
For the case of ℎ = 6.35 mm, the solution derived with 24 × 24 four-node shell elements by Sze et al. [31] is used as the major reference.The load-deflection curve is shown in Figure 18.It is seen that the solution in this study is very close to the solution of Sze et al. [31] in Figures 17 and 18.
In the third case of the present analysis, even a further reduction of the panel thickness was presumed taking ℎ = 3.175 mm.The load-deflection curve of very thin panel was presented in Figure 19.One can remark that the decrease of the thickness of the shell has reduced its stiffness; nevertheless, the range of rotations has not exceeded the limits of small rotations.the first is dominated by bending stiffness with large displacements; the second, at load level of  = 20, 000, is characterized by a very stiff response of the shell.The present converged solutions are compared with the results obtained by Sze et al. [31], who used the commercial code ABAQUS.The present results agree well with those obtained by Sze et al. [31] for the three curves considered.The deformed shape of the cylindrical shell, at load level of  = 40, 000, is presented in Figure 22.

Nonlinear Analysis of Pinched Cylindrical Shell.
Here, the cylindrical shell which is clamped at one end and subjected to a pair of concentrated loads at the other free end is studied.This example provides a severe test for finite element formulations.This example was considered in Brank et al. [29], Park et al. [36], Sansour and Kollmann [30], and Arciniega and Reddy [35] among others.
The two forces act in opposite directions as shown in Figure 23.Material properties are Young's modulus,  = 2.0685 × 10 7 , Poisson ratio, V = 0.3.The cylinder length is  = 304.8,and the radius  = 101.6 with thickness ℎ = 3.0.Onequarter of the structure is analyzed with 20×20 element mesh.In Figure 24, the radial displacements at points A and B of the shell are shown and compared with the reference solution of Sze et al. [31].Figure 24 shows that the load-deflection curve goes beyond the highest physically possible deflection of the loading point.The main purpose of the presented example was to test the behavior of the numerical model, although the authors are aware that the result of a real structure can differ from the presented solution due to a possible contact between two deformed regions of the shell in the postbuckling region.
The deformed shape of the semicylindrical shell subjected to an end pinching load, at load level of 0.35, is presented in Figure 25.

Concluding Remarks
In this study we have improved an 8-node ANS finite element model for linear and non-linear analysis of plate and shell structures.The accuracy and robustness of the improved formulation and its implementation are illustrated via a variety of shell geometries and deformations of shells.The objective of this paper is to present some results using the 8-node shell element when the sampling points for the strain components are changed.The 8-node ANS shell element using the sampling points earlier proposed (Bathe and Dvorkin [3], Lakshminarayana and Kailash [2], Bucalem and Bathe [4], Kim et al. [6], and Han et al., [8]) is considerably a powerful element and has shown good convergence behavior in many problem solutions.However, the research paper of predicted strain components show that the proper combination of sampling points is not predicted with excellent accuracy or reveals persistence of locking problems.
In order to improve the 8-node ANS shell element, a new combination of sampling points is adopted for linear analysis of isotropic and laminated composite structures.The present assumed strain method completely removes both membrane and shear locking behavior even when full integration is        used in the formulation.A new combination of the sampling points for the in-plane normal, in-plane shear, and transverse shear strain components results in significantly better results.Optimality in the convergence behavior is retained, and all strain components are predicted with reasonable accuracy.Other combinations of sampling points are also considered but do not result into a further improvement in the predictive capability of the element.From several numerical examples, the new combination 8-node ANS shell element shows better performance compared with other combination shell elements.The present combination of sampling points ( * 6 ) could be easily implemented into finite element code and used for the practical purpose.In particular, numerical   examples for annular plate presented herein clearly show the validity of the present approach and the accuracy of the developed nonlinear shell element.Future work will be useful to extend this work to dynamic and buckling analysis of isotropic and laminated composite shell structures.

Figure 3 :
Figure 3: Sampling points of the 8-node ANS shell element.
(a) For the Case of / = 1.0.

Figure 5 :
Figure 5: Boundary displacement conditions for patch tests.

Figure 8 :
Figure 8: Distorted mesh of simply supported and clamped rectangular plate.

Figure 14 :Figure 15 :
Figure 14: Load-deflection curves for the slit annular plate lifted by line force .

Figure 25 :
Figure 25: The deformed shape of pinched cylindrical shell at 0.35.

Table 1
[23,24]al tests of the present shell elements are carried out and validated using FEAP program (Zienkiewicz and Taylor[23,24]).All the results of the shell element showed very good agreement with references.Several examples are given to demonstrate the efficiency and accuracy of the present shell element.

Table 1 :
List of shell elements used for comparison.

Table 2 :
Results of patch test under bending (reference solution:
other patterns,  and  6 , still have the errors.Therefore, in this study, the improved patterns are used and the case of  * 6 shows more correct results.4.2.2.Distorted Mesh(/ = 1.0).In order to test the general applicability of the present formulation for bending, a thin Mathematical Problems in Engineering (c) In-plane tension test

Table 9 :
Results of simply supported plate with rectangular and distorted meshes (reference solutions:  3 = 4.062 and 11.60).