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

Authors: Woo-Young Jung; Sung-Cheon Han
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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

Figure1 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ξ3 (Rikards et al, [14]). Further well-known geometric relations for the shell in normal coordinates are also used.Let us assume that vectorP- characterizes the position of an arbitrary point of the shell in the initial reference state (see point B in Figure 2) and vector Q- is the position of the same point (point B′) in the deformed state. The position of a point at the mid-surface of the shell (point A) in the initial state is characterized by vector P, and position of the same point in the deformed state (point A′) is characterized by vector Q. Normal curvilinear coordinates ξi=[ξα,ξ3] at the mid-surface of the shell in the initial state are defined by the right-handed triad of the base vectors [aα,a3]. Here, a unit vector a3 is normal to the mid-surface of shell. Therefore (see Figure 2),
(1)P-(ξi)=P(ξα)+ξ3a3.Figure 2
Kinematics of the first-order shear deformation theory.Similarly, curvilinear coordinatesξ-i=[ξ-α,ξ-3] in the deformed state are defined by the triad of vectors [Aα,A3]. In the deformed state, the vector A3 may not be perpendicular to the mid-surface of shell.Vector functionQ- can be expanded in a Taylor’s series with respect to coordinate ξ3 normal to the mid-surface of shell as follows:
(2)Q-(ξi)=Q(ξα)+ξ3∇⊗Q·a3+L.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-surfaceu and displacement of an arbitrary point of shell u- are introduced. The following expressions can be written (see Figure 2) as
(3)Q-=P-+u-,Q=P+u.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
(4)u-(ξα)=u(ξα)+ξ3φ(ξα),
where φ is vector of rotation at the mid-surface of shell
(5)φ(ξα)=A3(ξα)-a3(ξα),
where
(6)A3=G·a3.The incremental form of the displacement field may be written in terms of the nodal incremental vectorΔua as
(7)Δu-(ξα)=∑a=18Na(ξα)[I3×3ξ3ha2Ha]Δu-a,
where Δu-a={Δu1a,Δu2a,Δu3a,Δϕ1a,Δϕ2a,Δϕ3a} and Ha is the rotational matrix of the normal vector at node a(8)Ha=haTRaΦaTA,
where TRa 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:(9)TRa=[cosϕ2cosϕ3-cosϕ2sinϕ3sinϕ2cosϕ1sinϕ3+sinϕ1sinϕ2cosϕ3cosϕ1cosϕ3-sinϕ1sinϕ2sinϕ3-sinϕ1cosϕ2sinϕ1sinϕ3-cosϕ1sinϕ2cosϕ3sinθ1cosϕ3+cosϕ1sinϕ2sinϕ3cosϕ1cosϕ2],Φa=[0a3a-a2a-a3a0a1aa2a-a1a0],TA=[cosϕ2cosϕ3sinϕ30-cosϕ2sinϕ3cosϕ30sinϕ201].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)ΔEm=BmΔu-,ΔEb=ξ3BbΔu-,ΔEs=BsΔu-.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 energyU of the shell represented as a three-dimensional body is given by the expression, where in curvilinear coordinates, the contravariant stress tensor Sij is contracted with the covariant strain tensor Eij(11)U=12∫VSijEijdV,wheredV is a volume element of the shell.The Hooke’s law for each layer can be written as(12)Sij=CijklEkl.Here,Cijkl is the component of the elasticity tensor in the mid-surface metrics. Substituting (12) into (11) as the strain energy U can be expressed as
(13)U=12∫A∫-h/2h/2EijTCijklEkldξ3dA.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:(14){RmbRs}=[Dmb00Ds]{EmbEs},where Rmb={Nξ1,Nξ2,Nξ1ξ2,Mξ1,Mξ2,Mξ1ξ2}, and Rs={Qξ1,Qξ2}.The membrane-bendingDmb and shear strain rigidity Ds have 6×6 and 2×2 matrices, respectively
(15)Dmb=∫-h/2h/2[C1ξ3C1ξ3C1ξ32C1]dξ3,
where C1 is the elastic constitutive coefficient
(16)Ds=ks∫-h/2h/2[G00G]dξ3,
where ks is the shear correction factor (ks=5/6) and G is the shear modulus. When the entire cross section is isotropic elastic, then the rigidity matrix reduces to the following form:
(17)De=[hC1000h312C1000hC2],
where C1 and C2 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.Sampling points of the 8-node ANS shell element.
(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])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 S1 and S2 (Bathe and Dvorkin [3]). The six cases of the combinations of various sampling points are used in the analysis.


## 2.1. Kinematics of Shell

Figure1 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ξ3 (Rikards et al, [14]). Further well-known geometric relations for the shell in normal coordinates are also used.Let us assume that vectorP- characterizes the position of an arbitrary point of the shell in the initial reference state (see point B in Figure 2) and vector Q- is the position of the same point (point B′) in the deformed state. The position of a point at the mid-surface of the shell (point A) in the initial state is characterized by vector P, and position of the same point in the deformed state (point A′) is characterized by vector Q. Normal curvilinear coordinates ξi=[ξα,ξ3] at the mid-surface of the shell in the initial state are defined by the right-handed triad of the base vectors [aα,a3]. Here, a unit vector a3 is normal to the mid-surface of shell. Therefore (see Figure 2),
(1)P-(ξi)=P(ξα)+ξ3a3.Figure 2
Kinematics of the first-order shear deformation theory.Similarly, curvilinear coordinatesξ-i=[ξ-α,ξ-3] in the deformed state are defined by the triad of vectors [Aα,A3]. In the deformed state, the vector A3 may not be perpendicular to the mid-surface of shell.Vector functionQ- can be expanded in a Taylor’s series with respect to coordinate ξ3 normal to the mid-surface of shell as follows:
(2)Q-(ξi)=Q(ξα)+ξ3∇⊗Q·a3+L.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-surfaceu and displacement of an arbitrary point of shell u- are introduced. The following expressions can be written (see Figure 2) as
(3)Q-=P-+u-,Q=P+u.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
(4)u-(ξα)=u(ξα)+ξ3φ(ξα),
where φ is vector of rotation at the mid-surface of shell
(5)φ(ξα)=A3(ξα)-a3(ξα),
where
(6)A3=G·a3.The incremental form of the displacement field may be written in terms of the nodal incremental vectorΔua as
(7)Δu-(ξα)=∑a=18Na(ξα)[I3×3ξ3ha2Ha]Δu-a,
where Δu-a={Δu1a,Δu2a,Δu3a,Δϕ1a,Δϕ2a,Δϕ3a} and Ha is the rotational matrix of the normal vector at node a(8)Ha=haTRaΦaTA,
where TRa 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:(9)TRa=[cosϕ2cosϕ3-cosϕ2sinϕ3sinϕ2cosϕ1sinϕ3+sinϕ1sinϕ2cosϕ3cosϕ1cosϕ3-sinϕ1sinϕ2sinϕ3-sinϕ1cosϕ2sinϕ1sinϕ3-cosϕ1sinϕ2cosϕ3sinθ1cosϕ3+cosϕ1sinϕ2sinϕ3cosϕ1cosϕ2],Φa=[0a3a-a2a-a3a0a1aa2a-a1a0],TA=[cosϕ2cosϕ3sinϕ30-cosϕ2sinϕ3cosϕ30sinϕ201].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)ΔEm=BmΔu-,ΔEb=ξ3BbΔu-,ΔEs=BsΔu-.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 energyU of the shell represented as a three-dimensional body is given by the expression, where in curvilinear coordinates, the contravariant stress tensor Sij is contracted with the covariant strain tensor Eij(11)U=12∫VSijEijdV,wheredV is a volume element of the shell.The Hooke’s law for each layer can be written as(12)Sij=CijklEkl.Here,Cijkl is the component of the elasticity tensor in the mid-surface metrics. Substituting (12) into (11) as the strain energy U can be expressed as
(13)U=12∫A∫-h/2h/2EijTCijklEkldξ3dA.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:(14){RmbRs}=[Dmb00Ds]{EmbEs},where Rmb={Nξ1,Nξ2,Nξ1ξ2,Mξ1,Mξ2,Mξ1ξ2}, and Rs={Qξ1,Qξ2}.The membrane-bendingDmb and shear strain rigidity Ds have 6×6 and 2×2 matrices, respectively
(15)Dmb=∫-h/2h/2[C1ξ3C1ξ3C1ξ32C1]dξ3,
where C1 is the elastic constitutive coefficient
(16)Ds=ks∫-h/2h/2[G00G]dξ3,
where ks is the shear correction factor (ks=5/6) and G is the shear modulus. When the entire cross section is isotropic elastic, then the rigidity matrix reduces to the following form:
(17)De=[hC1000h312C1000hC2],
where C1 and C2 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.Sampling points of the 8-node ANS shell element.
(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])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 S1 and S2 (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(18)∫δ(ΔLE)TCΔLEdV+∫Sδ(ΔNLE)dV=δt+ΔtWext-∫δ(ΔLE)TSdV,where superscriptt which is generally used as the current configuration is ignored in (18) and superscript t+Δt is the adjust incremented configuration; δt+ΔtWext is the external virtual work in t+Δt.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 (KL) is obtained as follows:
(19)∫δ(ΔLE)TCΔLΔEdV=δΔUT(∫BTCBdV)ΔU=δΔUTKLΔU.Finally, the element stiffness matrix has6×6 size on the reference surface of shell element
(20)[KL]=∫[KL11KL12KL21KL22]6×6dA,
where the submatrix of [KL] 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 ([KG]) is obtained:
(21)∫Sαβδ(ΔNLE)dV=∫δ(ΔNLE)TSαβdV=∫δΔΩTΔQTSαβdV.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 with6×6 submatrix as follows:
(22)[KG]=∫[KG11KG12KG21KG22]6×6dA,
where, the sub-matrix of [KG] is shown in Han et al. [10]. Then the final assembled incremental non-linear equilibrium equation can be written as
(23)([KL]+[KG])ΔU=F-t+Δt-F,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].


### 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 has6×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
(24)Ut=ktG×∫Ve[∂w1∂z2αt(ξ1,ξ2)-12{∂w2∂z1(ξ1,ξ2,0)-∂w1∂z2(ξ1,ξ2,0)}]2dV,
where kt is a parameter to be determined (the value of 0.1 suggested); G is the shear modulus; Ve is the volume of the element; and dV is the volume element. After integration throughout the thickness, (24) can be written as
(25)Ut=ktGh×∫Se[αt(ξ1,ξ2)∂w2∂z1-12{∂w2∂z1(ξ1,ξ2,0)-∂w1∂z2(ξ1,ξ2,0)}]2dξ1dξ2.IfktGh is chosen to be large relative to the factor Eh3 (which appears in the bending energy), (25) will play the role of penalty function and result in the desired constraint:
(26)αt(ξ1,ξ2)≈12{∂w2∂z1(ξ1,ξ2,0)-∂w1∂z2(ξ1,ξ2,0)}.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])
(27)Ut=u-TKtLu-.The torsional stiffness termKtL is added as described by Kanok-Nukulchai [22].


## 3.1. Linear Element Stiffness Matrix

If the strain displacement (8) is substituted into (18), the linearized element material stiffness matrix (KL) is obtained as follows:
(19)∫δ(ΔLE)TCΔLΔEdV=δΔUT(∫BTCBdV)ΔU=δΔUTKLΔU.Finally, the element stiffness matrix has6×6 size on the reference surface of shell element
(20)[KL]=∫[KL11KL12KL21KL22]6×6dA,
where the submatrix of [KL] 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 ([KG]) is obtained:
(21)∫Sαβδ(ΔNLE)dV=∫δ(ΔNLE)TSαβdV=∫δΔΩTΔQTSαβdV.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 with6×6 submatrix as follows:
(22)[KG]=∫[KG11KG12KG21KG22]6×6dA,
where, the sub-matrix of [KG] is shown in Han et al. [10]. Then the final assembled incremental non-linear equilibrium equation can be written as
(23)([KL]+[KG])ΔU=F-t+Δt-F,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].


## 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 has6×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
(24)Ut=ktG×∫Ve[∂w1∂z2αt(ξ1,ξ2)-12{∂w2∂z1(ξ1,ξ2,0)-∂w1∂z2(ξ1,ξ2,0)}]2dV,
where kt is a parameter to be determined (the value of 0.1 suggested); G is the shear modulus; Ve is the volume of the element; and dV is the volume element. After integration throughout the thickness, (24) can be written as
(25)Ut=ktGh×∫Se[αt(ξ1,ξ2)∂w2∂z1-12{∂w2∂z1(ξ1,ξ2,0)-∂w1∂z2(ξ1,ξ2,0)}]2dξ1dξ2.IfktGh is chosen to be large relative to the factor Eh3 (which appears in the bending energy), (25) will play the role of penalty function and result in the desired constraint:
(26)αt(ξ1,ξ2)≈12{∂w2∂z1(ξ1,ξ2,0)-∂w1∂z2(ξ1,ξ2,0)}.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])
(27)Ut=u-TKtLu-.The torsional stiffness termKtL 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 Table1. 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.
Name
Description (Patterns referred to Figure2)
Present (αδβ)
Sampling points for membrane (Patternα), in-plane shear (Pattern δ), and transverse shear (Pattern β)
Present (αδγ)
Sampling points for membrane (Patternα), in-plane shear (Pattern δ), and transverse shear (Pattern γ)
Present (βδγ)
Sampling points for membrane (Patternβ), in-plane shear (Pattern δ), and transverse shear (Pattern γ)
Present (βδγ*)
Sampling points for membrane (Patternβ), in-plane shear (Pattern δ), and transverse shear (Pattern γ*)
Present (βδγ6)
Sampling points for membrane (Patternβ), in-plane shear (Pattern δ), and transverse shear (Pattern γ6)
Present (βδγ6*)
Sampling points for membrane (Patternβ), in-plane shear (Pattern δ), and transverse shear (Pattern γ6*)
QUAD8
8-node shell element (MacNeal and Harder [7])
QUAD8*
8-node shell element (STRAND 7 [15])
QUAD8**
8-node ANS shell element (Lakshminarayana and Kailash [2])
XSHELL41
4-node quasiconforming shell element (Kim et al. [16])
XSHELL-8-ANS
8-node assumed natural strain shell element (XFINAS [17])
MITC8
8-node shell element using a mixed interpolation of tensorial components (Bathe and Dvorkin [3])


### 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:
(28)normalizedsolution=presentsolutionreferencesolution.Table 2
Results of patch test under bending (reference solution:θy=ML/EI=0.12×10-4).
Combinations of various sampling points
α
δ
β
α
δ
γ
β
δ
γ
β
δ
γ
*
β
δ
γ
6
β
δ
γ
6
*
Normalized solutions
1.000
1.000
1.000
1.000
1.000
1.000Table 3
Results of patch test under out-of-plane shear (reference solution:w=6SL/5GA=0.312×10-5).
Combinations of various sampling points
α
δ
β
α
δ
γ
β
δ
γ
β
δ
γ
*
β
δ
γ
6
β
δ
γ
6
*
Normalized solutions
0.053
0.989
0.989
0.989
0.989
1.000Table 4
Results of patch test under in-plane tension (reference solution:u=TL/EA=1.0×10-6).
Combinations of various sampling points
α
δ
β
α
δ
γ
β
δ
γ
β
δ
γ
*
β
δ
γ
6
β
δ
γ
6
*
Normalized solutions
1.000
0.965
0.992
0.992
0.992
0.992Figure 4
Mesh for patch test (Simo et al., [13]). Length of the square L=10; Young’s modulus E=1.0×107; Poisson’s ratio v=0.3; and thickness h=1.0.Boundary displacement conditions for patch tests.
(a)
Bending test
(b)
Transverse shear test
(c)
In-plane tension testLoading conditions for patch tests.
(a)
Bending test
(b)
Transverse shear test
(c)
In-plane tension testPatch test results indicate that various sampling points for assumed natural strain method in Tables2 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.


### 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 ofb/a=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 (b/a=1.0; reference solution: u3=4.062).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
6
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
1.009
1.000
1.000
1.000
1.016Table 6
Results of rectangular plate (b/a=1.0; reference solution: u3=5.60).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
6
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
1.057
1.001
1.001
0.997
1.090Table 7
Results of rectangular plate (b/a=5.0; reference solution: u3=12.97).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
6
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
1.009
1.000
1.000
1.000
1.010Table 8
Results of rectangular plate (b/a=5.0; reference solution: u3=7.23).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
6
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
0.845
0.845
0.897
0.975
0.867
8
×
4
1.029
0.991
0.991
—
—Table 9
Results of simply supported plate with rectangular and distorted meshes (reference solutions:u3=4.062 and 11.60).
Mesh types
Element
Uniform load(u3=4.062)
Concentrated load(u3=11.60)
Rectangular mesh
β
δ
γ
1.009
1.022
β
δ
γ
6
*
1.000
0.999
Distorted mesh
β
δ
γ
1.010
1.016
β
δ
γ
6
*
1.000
0.999Table 10
Results of clamped plate with rectangular and distorted meshes (reference solutions:u3=1.2637 and 5.60).
Mesh types
Elements
Uniform load(u3=1.2637)
Concentrated load(u3=5.60)
Rectangular mesh
Choi et al*
0.965
0.990
XSHELL41
0.997
0.992
β
δ
γ
1.052
1.057
β
δ
γ
6
*
1.002
1.001
Distorted mesh
Choi et al.*
0.976
1.011
XSHELL41
1.023
0.995
β
δ
γ
1.049
1.039
β
δ
γ
6
*
0.999
0.998
*Results of Choi et al. [20] and XSHELL41 are referenced from Kim et al. [16].


#### 4.2.1. Rectangular Mesh

(a) For the Case of  b/a=1.0.The reference vertical deflection at the center of the simply supported plate (shown in Figure7) 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  b/a=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 Tables7 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 γ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(b/a=1.0)

In order to test the general applicability of the present formulation for bending, a thin square plate,L/h=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 Tables9 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 βδγ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.


### 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 Figure9. 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 (h=0.1) under in-plane shear (reference solution: u2=0.08854).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8*
QUAD8**
6
×
1
0.991
0.991
0.991
0.869
0.942Table 12
Results of curved beam (h=0.1) under out-of-plane shear (reference solution: u3=0.5022).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8*
QUAD8**
6
×
1
0.970
0.972
0.972
0.958
0.962Figure 9
Curved Cantilever Beam.Results presented in Tables11 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 Figure10), 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 (h=0.32) under in-plane shear (reference solution: u2=5.424×10-3).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8
QUAD8**
6
×
1
0.997
0.997
0.993
0.998(12×2)
0.998(12×2)Table 14
Results of twisted beam (h=0.32) under out-of-plane shear (reference solution: u1=1.754×10-3).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8
QUAD8**
6
×
1
1.001
1.004
0.997
0.998(12×2)
0.998(12×2)Figure 10
Twisted cantilever beam.Results are presented in Tables13 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 Figure11. 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 thex-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:uA=0.093).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
0.542
0.543
0.516
0.832
0.387
6
×
6
0.915
0.915
0.902
1.003
0.773
8
×
8
0.984
0.984
0.978
C*
0.950
*The grade C means that the QUAD8 gives accuracy 20% > error > 10% (MacNeal and Harder [7]).This indicates that the results obtained by our elements are more accurate than the solution by the QUAD8** element. Note that the error of8×8 mesh for the QUAD8 element is greater than 10%.


### 4.6. Pinched Cylinder with End Diaphragm

This problem, shown in Figure12, 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.Table 16
Results of pinched cylinder with end diaphragm (reference solution:uA=1.8248×10-5).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
XSHELL-8-ANS
MITC8
4
×
4
0.687
0.720
0.781
0.738
0.833 (3×3)
8
×
8
1.006
1.024
0.964
0.945
0.990Figure 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 βδγ6* 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 Figure13). 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 forcesF.The line forceF 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 F=1.45.Figure 14
Load-deflection curves for the slit annular plate lifted by line forceF.Figure 15
The deformed shape of slip annular plate atF=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 with12×12 mesh sizes, as shown in Figure 16. Three different thickness values are considered; that is, h=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 L=508 mm, curvature of the circular edges (1/R) is 1/2,540 mm, and the circumferential length is 508 mm. The material properties are Young’s modulus E=3,102.75N/mm2, and Poisson’s ratio v=0.3.Figure 16
Hinged cylindrical shell; curved edges are free, and straight edges are hinged(P=1,000).For the case ofh= 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 (h=12.7mm).For the case ofh= 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 (h=6.35mm).In the third case of the present analysis, even a further reduction of the panel thickness was presumed takingh=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 (h=3.175mm).


### 4.9. Nonlinear Analysis of Open-Ended Cylindrical Shell

The large deformation analysis of a stretched cylinder, as shown in Figure20, 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 E=10.5×106 and v=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.(P=40,000).The radial displacements at points A, B, C, and D of the shell are shown in Figure21. 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 P=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 P=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 atP=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 Figure23. Material properties are Young’s modulus, E=2.0685×107, Poisson ratio, v=0.3. The cylinder length is L=304.8, and the radius R=101.6 with thickness h=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.35P, is presented in Figure 25.Figure 23
Problem definition and an initial configuration(P=2,000).Figure 24
Load-deflection curves for the pinched cylindrical shell.Figure 25
The deformed shape of pinched cylindrical shell at 0.35P.


## 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:
(28)normalizedsolution=presentsolutionreferencesolution.Table 2
Results of patch test under bending (reference solution:θy=ML/EI=0.12×10-4).
Combinations of various sampling points
α
δ
β
α
δ
γ
β
δ
γ
β
δ
γ
*
β
δ
γ
6
β
δ
γ
6
*
Normalized solutions
1.000
1.000
1.000
1.000
1.000
1.000Table 3
Results of patch test under out-of-plane shear (reference solution:w=6SL/5GA=0.312×10-5).
Combinations of various sampling points
α
δ
β
α
δ
γ
β
δ
γ
β
δ
γ
*
β
δ
γ
6
β
δ
γ
6
*
Normalized solutions
0.053
0.989
0.989
0.989
0.989
1.000Table 4
Results of patch test under in-plane tension (reference solution:u=TL/EA=1.0×10-6).
Combinations of various sampling points
α
δ
β
α
δ
γ
β
δ
γ
β
δ
γ
*
β
δ
γ
6
β
δ
γ
6
*
Normalized solutions
1.000
0.965
0.992
0.992
0.992
0.992Figure 4
Mesh for patch test (Simo et al., [13]). Length of the square L=10; Young’s modulus E=1.0×107; Poisson’s ratio v=0.3; and thickness h=1.0.Boundary displacement conditions for patch tests.
(a)
Bending test
(b)
Transverse shear test
(c)
In-plane tension testLoading conditions for patch tests.
(a)
Bending test
(b)
Transverse shear test
(c)
In-plane tension testPatch test results indicate that various sampling points for assumed natural strain method in Tables2 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.


## 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 ofb/a=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 (b/a=1.0; reference solution: u3=4.062).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
6
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
1.009
1.000
1.000
1.000
1.016Table 6
Results of rectangular plate (b/a=1.0; reference solution: u3=5.60).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
6
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
1.057
1.001
1.001
0.997
1.090Table 7
Results of rectangular plate (b/a=5.0; reference solution: u3=12.97).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
6
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
1.009
1.000
1.000
1.000
1.010Table 8
Results of rectangular plate (b/a=5.0; reference solution: u3=7.23).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
6
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
0.845
0.845
0.897
0.975
0.867
8
×
4
1.029
0.991
0.991
—
—Table 9
Results of simply supported plate with rectangular and distorted meshes (reference solutions:u3=4.062 and 11.60).
Mesh types
Element
Uniform load(u3=4.062)
Concentrated load(u3=11.60)
Rectangular mesh
β
δ
γ
1.009
1.022
β
δ
γ
6
*
1.000
0.999
Distorted mesh
β
δ
γ
1.010
1.016
β
δ
γ
6
*
1.000
0.999Table 10
Results of clamped plate with rectangular and distorted meshes (reference solutions:u3=1.2637 and 5.60).
Mesh types
Elements
Uniform load(u3=1.2637)
Concentrated load(u3=5.60)
Rectangular mesh
Choi et al*
0.965
0.990
XSHELL41
0.997
0.992
β
δ
γ
1.052
1.057
β
δ
γ
6
*
1.002
1.001
Distorted mesh
Choi et al.*
0.976
1.011
XSHELL41
1.023
0.995
β
δ
γ
1.049
1.039
β
δ
γ
6
*
0.999
0.998
*Results of Choi et al. [20] and XSHELL41 are referenced from Kim et al. [16].


### 4.2.1. Rectangular Mesh

(a) For the Case of  b/a=1.0.The reference vertical deflection at the center of the simply supported plate (shown in Figure7) 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  b/a=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 Tables7 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 γ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(b/a=1.0)

In order to test the general applicability of the present formulation for bending, a thin square plate,L/h=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 Tables9 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 βδγ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.


## 4.2.1. Rectangular Mesh

(a) For the Case of  b/a=1.0.The reference vertical deflection at the center of the simply supported plate (shown in Figure7) 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  b/a=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 Tables7 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 γ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(b/a=1.0)

In order to test the general applicability of the present formulation for bending, a thin square plate,L/h=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 Tables9 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 βδγ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.


## 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 Figure9. 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 (h=0.1) under in-plane shear (reference solution: u2=0.08854).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8*
QUAD8**
6
×
1
0.991
0.991
0.991
0.869
0.942Table 12
Results of curved beam (h=0.1) under out-of-plane shear (reference solution: u3=0.5022).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8*
QUAD8**
6
×
1
0.970
0.972
0.972
0.958
0.962Figure 9
Curved Cantilever Beam.Results presented in Tables11 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 Figure10), 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 (h=0.32) under in-plane shear (reference solution: u2=5.424×10-3).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8
QUAD8**
6
×
1
0.997
0.997
0.993
0.998(12×2)
0.998(12×2)Table 14
Results of twisted beam (h=0.32) under out-of-plane shear (reference solution: u1=1.754×10-3).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8
QUAD8**
6
×
1
1.001
1.004
0.997
0.998(12×2)
0.998(12×2)Figure 10
Twisted cantilever beam.Results are presented in Tables13 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 Figure11. 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 thex-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:uA=0.093).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
QUAD8
QUAD8**
4
×
4
0.542
0.543
0.516
0.832
0.387
6
×
6
0.915
0.915
0.902
1.003
0.773
8
×
8
0.984
0.984
0.978
C*
0.950
*The grade C means that the QUAD8 gives accuracy 20% > error > 10% (MacNeal and Harder [7]).This indicates that the results obtained by our elements are more accurate than the solution by the QUAD8** element. Note that the error of8×8 mesh for the QUAD8 element is greater than 10%.


## 4.6. Pinched Cylinder with End Diaphragm

This problem, shown in Figure12, 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.Table 16
Results of pinched cylinder with end diaphragm (reference solution:uA=1.8248×10-5).
Element size
Normalized solutions
β
δ
γ
β
δ
γ
*
β
δ
γ
6
*
XSHELL-8-ANS
MITC8
4
×
4
0.687
0.720
0.781
0.738
0.833 (3×3)
8
×
8
1.006
1.024
0.964
0.945
0.990Figure 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 βδγ6* 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 Figure13). 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 forcesF.The line forceF 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 F=1.45.Figure 14
Load-deflection curves for the slit annular plate lifted by line forceF.Figure 15
The deformed shape of slip annular plate atF=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 with12×12 mesh sizes, as shown in Figure 16. Three different thickness values are considered; that is, h=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 L=508 mm, curvature of the circular edges (1/R) is 1/2,540 mm, and the circumferential length is 508 mm. The material properties are Young’s modulus E=3,102.75N/mm2, and Poisson’s ratio v=0.3.Figure 16
Hinged cylindrical shell; curved edges are free, and straight edges are hinged(P=1,000).For the case ofh= 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 (h=12.7mm).For the case ofh= 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 (h=6.35mm).In the third case of the present analysis, even a further reduction of the panel thickness was presumed takingh=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 (h=3.175mm).


## 4.9. Nonlinear Analysis of Open-Ended Cylindrical Shell

The large deformation analysis of a stretched cylinder, as shown in Figure20, 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 E=10.5×106 and v=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.(P=40,000).The radial displacements at points A, B, C, and D of the shell are shown in Figure21. 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 P=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 P=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 atP=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 Figure23. Material properties are Young’s modulus, E=2.0685×107, Poisson ratio, v=0.3. The cylinder length is L=304.8, and the radius R=101.6 with thickness h=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.35P, is presented in Figure 25.Figure 23
Problem definition and an initial configuration(P=2,000).Figure 24
Load-deflection curves for the pinched cylindrical shell.Figure 25
The deformed shape of pinched cylindrical shell at 0.35P.


## 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 (βδγ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.



---

*Source: 276304-2013-08-06.xml*