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

Abstract

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.

1. 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 thin-walled structures, the accuracy of the element showed less than the predicted results.

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.

2. 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.

Figure 1 

Geometry of 8-node 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 𝜉₃ (Rikards et al, [14]). Further well-known geometric relations for the shell in normal coordinates are also used.

Let us assume that vector 𝐏 characterizes the position of an arbitrary point of the shell in the initial reference state (see point 𝐵 in Figure 2) and vector 𝐐 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 𝐏, and position of the same point in the deformed state (point 𝐴^′) is characterized by vector 𝐐. Normal curvilinear coordinates 𝜉^𝑖=[𝜉^𝛼,𝜉₃] at the mid-surface of the shell in the initial state are defined by the right-handed triad of the base vectors [𝐚_𝛼,𝐚₃]. Here, a unit vector 𝐚₃ is normal to the mid-surface of shell. Therefore (see Figure 2), 𝐏𝜉^𝑖=𝐏𝜉^𝛼+𝜉₃𝐚₃. (1)

Figure 2 

Kinematics of the first-order shear deformation theory.

Similarly, curvilinear coordinates 𝜉^𝑖=[𝜉^𝛼,𝜉₃] in the deformed state are defined by the triad of vectors [𝐀_𝛼,𝐀₃]. In the deformed state, the vector 𝐀₃ may not be perpendicular to the mid-surface of shell.

Vector function 𝐐 can be expanded in a Taylor’s series with respect to coordinate 𝜉₃ normal to the mid-surface of shell as follows: 𝐐𝜉^𝑖=𝐐𝜉^𝛼+𝜉₃∇⊗𝐐⋅𝐚₃+L. (2)

Here, ∇⊗𝐐=𝐀_𝛼⊗𝐚^𝛼=𝐆 is a second-order tensor (see Luri’e, [18]), which characterizes the gradient of strains with respect to coordinate 𝜉₃, where ∇ is a Hamiltonian operator and ⊗ is a dyadic product of the tensors.

Further, the notations for displacement at the mid-surface 𝐮 and displacement of an arbitrary point of shell 𝐮 are introduced. The following expressions can be written (see Figure 2) as 𝐐=𝐏+𝐮,𝐐=𝐏+𝐮. (3)

From expressions (1)–(3), the representation of the displacement 𝐮 of an arbitrary point of the shell for the first-order approximation can be obtained from 𝐮𝜉^𝛼=𝐮𝜉^𝛼+𝜉₃𝝋𝜉^𝛼, (4) where 𝝋 is vector of rotation at the mid-surface of shell 𝝋𝜉^𝛼=𝐀₃𝜉^𝛼−𝐚₃𝜉^𝛼, (5) where 𝐀₃=𝐆⋅𝐚₃. (6)

The incremental form of the displacement field may be written in terms of the nodal incremental vector Δ𝐮^𝑎 as Δ𝐮𝜉^𝛼=⁸_𝑎=1𝑁^𝑎𝜉^𝛼𝐈_3×3𝜉₃ℎ^𝑎2𝐇^𝑎Δ𝐮^𝑎, (7) where Δ𝑢^𝑎={Δ𝑢^𝑎₁,Δ𝑢^𝑎₂,Δ𝑢^𝑎₃,Δ𝜙^𝑎₁,Δ𝜙^𝑎₂,Δ𝜙^𝑎₃} and 𝐇^𝑎 is the rotational matrix of the normal vector at node 𝑎𝐇^𝑎=ℎ^𝑎𝐓^𝑎_𝐑Φ^𝑎𝐓_𝐀, (8) where 𝐓^𝑎_𝐑 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 𝜙₁, 𝜙₂, and 𝜙₃ as follows:𝐓^𝑎_𝐑=cos𝜙₂cos𝜙₃−cos𝜙₂sin𝜙₃sin𝜙₂cos𝜙₁sin𝜙₃+sin𝜙₁sin𝜙₂cos𝜙₃cos𝜙₁cos𝜙₃−sin𝜙₁sin𝜙₂sin𝜙₃−sin𝜙₁cos𝜙₂sin𝜙₁sin𝜙₃−cos𝜙₁sin𝜙₂cos𝜙₃sin𝜃₁cos𝜙₃+cos𝜙₁sin𝜙₂sin𝜙₃cos𝜙₁cos𝜙₂,Φ^𝑎=0𝐚^𝑎₃−𝐚^𝑎₂−𝐚^𝑎₃0𝐚^𝑎₁𝐚^𝑎₂−𝐚^𝑎₁0,𝐓_𝐀=cos𝜙₂cos𝜙₃sin𝜙₃0−cos𝜙₂sin𝜙₃cos𝜙₃0sin𝜙₂.01 (9)

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: Δ𝐸^𝑚=𝐁_𝐦Δ𝐮,Δ𝐸^𝑏=𝜉₃𝐁_𝐛Δ𝐮,Δ𝐸^𝑠=𝐁_𝐬Δ𝐮. (10)

The shell theory outlined above is the so-called first-order shear deformation theory with six degrees of freedom.

2.2. 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 𝐸_𝑖𝑗1𝑈=2_𝑉𝑆^𝑖𝑗𝐸_𝑖𝑗𝑑𝑉, (11)

where 𝑑𝑉 is a volume element of the shell.

The Hooke’s law for each layer can be written as 𝑆^𝑖𝑗=𝐶^{𝑖𝑗𝑘𝑙}𝐸_𝑘𝑙. (12)

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 1𝑈=2_𝐴^ℎ/2_−ℎ/2𝐸^𝑇_𝑖𝑗𝐶^{𝑖𝑗𝑘𝑙}𝐸_𝑘𝑙𝑑𝜉₃𝑑𝐴. (13)

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: 𝐑_𝑚𝑏𝐑_𝑠=𝐃_𝑚𝑏00𝐃_𝑠𝐸^𝑚𝑏𝐸^𝑠, (14)where 𝐑_𝑚𝑏={𝑁_𝜉¹,𝑁_𝜉²,𝑁_𝜉¹_𝜉²,𝑀_𝜉¹,𝑀_𝜉²,𝑀_𝜉¹_𝜉²}, and 𝐑_𝑠={𝑄_𝜉¹,𝑄_𝜉²}.

The membrane-bending 𝐃_𝑚𝑏 and shear strain rigidity 𝐃_𝑠 have 6×6 and 2×2 matrices, respectively 𝐃_𝑚𝑏=^ℎ/2_−ℎ/2𝐂_𝟏𝜉₃𝐂_𝟏𝜉₃𝐂_𝟏𝜉²₃𝐂_𝟏𝑑𝜉₃, (15) where 𝐂₁ is the elastic constitutive coefficient 𝐃_𝑠=𝑘_𝑠^ℎ/2_−ℎ/2𝐺00𝐺𝑑𝜉₃, (16) 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: 𝐃_𝑒=ℎ𝐂₁0ℎ00₃𝐂12¹000ℎ𝐂₂, (17) where 𝐂₁ and 𝐂₂ are the elastic constitutive coefficients of membrane-bending and transverse shear parts.

2.3. Various Enhanced Strain Interpolation

Many researchers addressed 8-node shell finite elements (Hinton and Huang [1]; Lakshminarayana and Kailash [2]; MacNeal and Harder [19]; Bucalem and Bathe [4]; Kim et al. [5, 6]; Han et al. [8]). 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.

(a) Pattern 𝛼 (following Bathe et al. [9])

(b) Pattern 𝛽 (following Han et al. [10])

(c) Pattern 𝛿 (following Han et al. [10])

(d) Patterns 𝛾,𝛾6 (following Hwang [11])

(e) Patterns 𝛾∗,𝛾∗6 (following Polit et al. [12])

Figure 3 

Sampling points of the 8-node ANS shell element.

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 𝛾^∗₆ patterns. In the 𝛾₆ and 𝛾^∗₆ patterns, the strain component of center point is replaced by the mean of the components at points 𝑆₁ and 𝑆₂ (Bathe and Dvorkin [3]). The six cases of the combinations of various sampling points are used in the analysis.

3. Incremental Equation of Equilibrium

The generalized Hook’s law at large strain does not represent an approximate material behavior description because stress-strain 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 𝛿Δ^𝐿𝐸^𝑇𝐂Δ^𝐿𝐸𝑑𝑉+𝐒𝛿Δ^𝑁𝐿𝐸=𝑑𝑉_𝑡+Δ𝑡𝛿𝑊_ext−𝛿Δ^𝐿𝐸^𝑇𝐒𝑑𝑉, (18)

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.

3.1. Linear Element Stiffness Matrix

If the strain displacement (8) is substituted into (18), the linearized element material stiffness matrix (𝐊_𝐋) is obtained as follows: 𝛿Δ^𝐿𝐸^𝑇𝐂Δ^𝐿Δ𝐸𝑑𝑉=𝛿Δ𝐔^𝐓𝐁^𝐓𝐂𝐁𝑑𝑉Δ𝐔=𝛿Δ𝐔^𝐓𝐊_𝐋Δ𝐔. (19)

Finally, the element stiffness matrix has 6×6 size on the reference surface of shell element 𝐊_𝐋𝐊=^𝐋𝟏𝟏𝐊_𝐋^𝟏𝟐𝐊_𝐋^𝟐𝟏𝐊_𝐋^𝟐𝟐_6×6𝑑𝐴, (20) where the submatrix of [𝐊_𝐋] is shown in Han et al. [10].

3.2. 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 ([𝐊_𝐆]) is obtained: 𝑆_𝛼𝛽𝛿Δ^𝑁𝐿𝐸𝑑𝑉=𝛿Δ^𝑁𝐿𝐸^𝑇𝑆_𝛼𝛽𝑑𝑉=𝛿Δ𝛀^𝐓Δ𝐐^𝐓𝑆_𝛼𝛽𝑑𝑉. (21)

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: 𝐊_𝐆𝐊=^𝐆𝟏𝟏𝐊_𝐆^𝟏𝟐𝐊_𝐆^𝟐𝟏𝐊_𝐆^𝟐𝟐_6×6𝑑𝐴, (22) where, the sub-matrix of [𝐊_𝐆] is shown in Han et al. [10]. Then the final assembled incremental non-linear equilibrium equation can be written as 𝐊_𝐋+𝐊_𝐆Δ𝐔=^𝑡+Δ𝑡𝐅−𝐅, (23)

where, 𝐅 and 𝐅 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].

3.3. 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 𝑈_𝑡=𝑘_𝑡𝐺×_𝑉𝑒𝛼_𝑡𝜉₁,𝜉₂−12𝜕𝑤₂𝜕𝑧₁𝜉₁,𝜉₂,0−𝜕𝑤₁𝜕𝑧₂𝜉₁,𝜉₂,0²𝑑𝑉, (24) 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 𝑈_𝑡=𝑘_𝑡×𝐺ℎ_𝑆𝑒𝛼_𝑡𝜉₁,𝜉₂−12𝜕𝑤₂𝜕𝑧₁𝜉₁,𝜉₂,0−𝜕𝑤₁𝜕𝑧₂𝜉₁,𝜉₂,0²𝑑𝜉₁𝑑𝜉₂. (25)

If 𝑘_𝑡𝐺ℎ is chosen to be large relative to the factor 𝐸ℎ³ (which appears in the bending energy), (25) will play the role of penalty function and result in the desired constraint: 𝛼_𝑡𝜉₁,𝜉₂≈12𝜕𝑤₂𝜕𝑧₁𝜉₁,𝜉₂,0−𝜕𝑤₁𝜕𝑧₂𝜉₁,𝜉₂,0. (26)

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, 𝐮, by shape functions. This gives (25) in the form (Kanok-Nukulchai [22]) 𝑈_𝑡=𝐮^𝐓𝐊_𝐭𝐋𝐮. (27)

The torsional stiffness term 𝐊_𝐭𝐋 is added as described by Kanok-Nukulchai [22].

4. Numerical Results of the 8-Node Shell Element

In this section, we describe numerical tests of improved 8-node shell elements listed in Table 1. They are based on the various combinations of assumed natural strain interpolation functions. Several numerical examples are solved to conduct the performance of the shell element in linear and non-linear applications. The patch test, distortion test, and other various numerical 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.

4.1. 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. The normalized solutions are presented in the nondimensional form using the following equation: normalizedsolution=presentsolution.referencesolution (28)

Table 2 

Results of patch test under bending (reference solution: 𝜃_𝑦=𝑀𝐿/𝐸𝐼=0.12×10⁻⁴).

Table 3 

Results of patch test under out-of-plane shear (reference solution: 𝑤=6𝑆𝐿/5𝐺𝐴=0.312×10⁻⁵).

Table 4 

Results of patch test under in-plane tension (reference solution: 𝑢=𝑇𝐿/𝐸𝐴=1.0×10⁻⁶).

Figure 4 

Mesh for patch test (Simo et al., [13]). Length of the square 𝐿=10; Young’s modulus 𝐸=1.0×107; Poisson’s ratio 𝑣=0.3; and thickness ℎ=1.0.

(a) Bending test  

(b) Transverse shear test

(c) In-plane tension test

Figure 5 

Boundary displacement conditions for patch tests.

(a) Bending test

(b) Transverse shear test

(c) In-plane tension test

Figure 6 

Loading conditions for patch tests.

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 𝛼𝛿𝛾, 𝛽𝛿𝛾, 𝛽𝛿𝛾₆, and 𝛽𝛿𝛾^∗₆ 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.

4.2. 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.

Table 5 

Results of rectangular plate (𝑏/𝑎=1.0; reference solution: 𝑢₃=4.062).

Table 6 

Results of rectangular plate (𝑏/𝑎=1.0; reference solution: 𝑢₃=5.60).

Table 7 

Results of rectangular plate (𝑏/𝑎=5.0; reference solution: 𝑢₃=12.97).

Table 8 

Results of rectangular plate (𝑏/𝑎=5.0; reference solution: 𝑢₃=7.23).

Table 9 

Results of simply supported plate with rectangular and distorted meshes (reference solutions: 𝑢₃=4.062 and 11.60).

Table 10 

Results of clamped plate with rectangular and distorted meshes (reference solutions: 𝑢₃=1.2637 and 5.60).

4.2.1. Rectangular Mesh

(a) For the Case of  𝑏/𝑎=1.0.

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.(1)Simply supported plate—uniform load.(2)Clamped plate—concentrated load.

Figure 7 

Simply supported and clamped rectangular plate.

(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.(1)Simply supported plate—uniform load.(2)Clamped plate—concentrated load.

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 for the application of 𝛽 pattern in transverse shear strain. A specifically designed application must be developed in order to remove these errors in the 8-node shell element because other patterns, 𝛾 and 𝛾₆, still have the errors. Therefore, in this study, the improved patterns are used and the case of 𝛽𝛿𝛾^∗₆ 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 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].

Figure 8 

Distorted mesh of simply supported and clamped rectangular plate.

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.(1)Simply supported plate—uniform and concentrated load.(2)Clamped plate—uniform and concentrated load.

The cases of 𝛽𝛿𝛾 and 𝛽𝛿𝛾^∗₆ 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.

4.3. 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-of-plane 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 Tables 11-12.

Table 11 

Results of curved beam (ℎ=0.1) under in-plane shear (reference solution: 𝑢₂=0.08854).

Table 12 

Results of curved beam (ℎ=0.1) under out-of-plane shear (reference solution: 𝑢₃=0.5022).

Figure 9 

Curved Cantilever Beam.

Results presented in Tables 11 and 12 show the outstanding performance of the proposed element.

4.4. 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]).

Table 13 

Results of twisted beam (ℎ=0.32) under in-plane shear (reference solution: 𝑢₂=5.424×10⁻³).

Table 14 

Results of twisted beam (ℎ=0.32) under out-of-plane shear (reference solution: 𝑢₁=1.754×10⁻³).

Figure 10 

Twisted cantilever beam.

Results are presented in Tables 13 and 14, showing more correct results than the result of references in spite of the 6×1 mesh.

4.5. 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.

Figure 11 

Definition of loading points and initial configuration of hemispherical shell.

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.

Table 15 

Results of hemispherical shell with 18° hole (reference solution: 𝑢_𝐴=0.093).

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%.

4.6. Pinched Cylinder with End Diaphragm

This problem, shown in Figure 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⁻⁴ (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.

Table 16 

Results of pinched cylinder with end diaphragm (reference solution: 𝑢_𝐴=1.8248×10⁻⁵).

Figure 12 

Problem of pinched cylinder with end diaphragms.

The cases of 𝛽𝛿𝛾 and 𝛽𝛿𝛾^∗ show that the results are more conservative than the reference solution. However, the solution in the case of 𝛽𝛿𝛾^∗₆ is more accurate than that in the reference element XSHELL-8-ANS.

4.7. 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.

Figure 13 

The slit annular plate loaded with various end forces 𝐹.

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 arc-length 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.

Figure 14 

Load-deflection curves for the slit annular plate lifted by line force 𝐹.

Figure 15 

The deformed shape of slip annular plate at 𝐹=1.45.

4.8. 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⁻⁶ 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.75N/mm², and Poisson’s ratio 𝑣=0.3.

Figure 16 

Hinged cylindrical shell; curved edges are free, and straight edges are hinged (𝑃=1,000).

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.

Figure 17 

Load-deflection curves for the hinged shell (ℎ=12.7mm).

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.

Figure 18 

Load-deflection curves for the hinged shell (ℎ=6.35mm).

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.

Figure 19 

Load-deflection curves for the hinged shell (ℎ=3.175mm).

4.9. Nonlinear Analysis of Open-Ended Cylindrical Shell

The large deformation analysis of a stretched cylinder, as shown in Figure 20, is carried out with free ends. It is subjected to a pair of concentrated forces. One octant of the cylinder is analyzed with 16×12 mesh sizes. Material properties are 𝐸=10.5×10⁶ and 𝑣=0.3125. The cylinder length is 10.35, the radius is 4.953, and the thickness is 0.094.

Figure 20 

Problem definition and an initial configuration. (𝑃=40,000).

The radial displacements at points A, B, C, and D of the shell are shown in Figure 21. It is noticed that the load-deflection curves of the shell have two different regions: 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.

Figure 21 

Load-deflection curves for the open-ended cylindrical shell.

Figure 22 

The deformed shape of open-ended shell at 𝑃=40,000.

4.10. 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⁷, Poisson ratio, 𝑣=0.3. The cylinder length is 𝐿=304.8, and the radius 𝑅=101.6 with thickness ℎ=3.0. One-quarter 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.

Figure 23 

Problem definition and an initial configuration (𝑃=2,000).

Figure 24 

Load-deflection curves for the pinched cylindrical shell.

Figure 25 

The deformed shape of pinched cylindrical shell at 0.35𝑃.

5. 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 (𝛽𝛿𝛾^∗₆) 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.

Acknowledgment

This work was supported by the National Research Foundation of Korea grant funded (NRF) by the Korean Government (MEST) (no. 2012-0008762).