A refined element-based Lagrangian shell element for geometrically nonlinear analysis of shell structures

Abstract

For the solution of geometrically nonlinear analysis of plates and shells, the formulation of a nonlinear nine-node refined first-order shear deformable element-based Lagrangian shell element is presented. Natural co-ordinate-based higher order transverse shear strains are used in present shell element. Using the 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, a refined first-order shear deformation theory for thin and thick shells, which results in parabolic through-thickness distribution of the transverse shear strains from the formulation based on the third-order shear deformation theory, is proposed. This formulation eliminates the need for shear correction factors in the first-order theory. To avoid difficulties resulting from large increments of the rotations, a scheme of attached reference system is used for the expression of rotations of shell normal. Numerical examples demonstrate that the present element behaves reasonably satisfactorily either for the linear or for geometrically nonlinear analysis of thin and thick plates and shells with large displacement but small strain. Especially, the nonlinear results of slit annular plates with various loads provided the benchmark to test the accuracy of related numerical solutions.

Introduction

Due to their attractive properties, shells are used in a variety of engineering structures, such as arch dams, dome structures, pressure-vessels, ship hulls, and aircraft fuel-ages. In order to fully utilize their special structural features, efficient and reliable analysis tools are needed for the prediction of their structural behaviors. Now the finite element method is a major tool for numerical analysis of shell structures because of their intrinsic complexities. Over the past 50 years, many shell elements have been developed. Many shell element formulations have been based on the so-called degenerate models introduced by Ahmad et al. ¹ While such elements are capable of dealing with thick plate and shell problems, their performance deteriorates rapidly as the element thickness becomes thin, which is called shear locking. In order to overcome the shear locking problems, Huang and Hinton ² developed a nine-node assumed strain shell element using the enhanced interpolation of the transverse shear strains in the natural coordinate system. The assumed natural strain method can be traced to MacNeal ³ and Dvorkin and Bathe. ⁴ Many finite elements employing the assumed strain method were then reported by other research groups independently.^5,6 Furthermore, some applications of the assumed strain method have been made in the field of large deformation analysis. In particular, assumed strain has been expressed in the form of natural strain because of its simplicity. In order to avoid shear locking in the plate and shell elements, Nguyen-Xuan et al. ⁷ and Nguyen-Thanh et al. ⁸ proposed a smoothed finite element, respectively. Recent years, Rodrigues et al. ⁹ developed smoothed finite element for analysis of composite plates. In this study, the element-based Lagrangian formulation (ELF) by Kanok-Nukulchai and Wong ¹⁰ is used and it makes implementation simpler and easier than the traditional Lagrangian formulations, especially when the assumed natural strain method is involved.

Zhang and Kim ¹¹ developed a refined non-conforming triangular plate/shell element for linear and geometrically nonlinear analysis of plates and shells. The drilling rotation constraint (RC) equation was used in the work ¹² with the Hu–Washizu (HW) functional. Various forms of governing functionals modified by the RC are presented in Wisniewski. ¹² Defining the assumed representations of stress and strain in the four-node HW elements, Wisniewski and Turska^13,14 used the skew coordinates instead of the natural ones. Firl and Bletzinger ¹⁵ studied the shape optimization of thin-walled structures. Shin and Lee ¹⁶ proposed a three-node triangular flat shell element based on the assumed natural deviatoric strain with drilling degrees of freedom. Also, a three-dimensional (3D) nonlinear analysis using a resultant eight-node solid element is presented by Kim et al. ¹⁷ and Ahmed et al. ¹⁸

The first-order shear deformable shell elements are usually capable of also producing good results for the in-plane strain and stress distribution, but their formulation results in constant transverse shear strains as opposed to the realistic parabolic distribution. As a result, the traction conditions at the shell surfaces are violated. They also require shear correction coefficients to correct the corresponding strain energy terms and these coefficients are problem dependent and are not always easy to determine. Numerous efforts have been made to overcome this disadvantage of the first-order formulation, most of which result in a higher order shear deformation theory.^19–21 The equilibrium approach suggested by Rolfes and Rohwer ²² was used for the calculation of the transverse shear stiffness in Kim at al. ⁵ Tanov and Tabiei ²³ presented an approach for treating the transverse shear strains and stresses in homogeneous shells resulting in parabolic distribution for both strains and stresses and, therefore, it eliminates the need of any shear correction factors. It requires only minor changes in the first-order shear deformation formulation. The implementation in shell element was shown to be quite simple and straightforward in Han et al.^24,25 and Jung and Han. ²⁶ The proposed nine-node shell element will consider parabolic through-thickness distribution of the refined transverse shear stresses and strains based on simple corrections to the first-order theory proposed by Han et al. ²⁴

In large deformation analysis, two steps are involved in obtaining the linearized nonlinear equations. First, nonlinear equations of the structural system have to be derived via Lagrangian formulations. Second, the linearization of the nonlinear equations is required to solve nonlinear equations by an iterative method such as the arc-length control method. So far, traditional Lagrangian formulations, which are total Lagrangian formulation (TLF) and updated Lagrangian formulation (ULF), have been extensively used for deriving the nonlinear equation. Unlike two formulations of the above, new Lagrangian formulation was introduced and it was referred to as the ELF since the fundamental element serves as a deformation reference in ELF. ⁷ It means that all equations governing a deformed body can be expressed with respect to natural coordinate systems. For the large deformation analysis of shells, the expression of rotations plays an important role to define a proper position of shell normal vector. In this study, three successive finite rotations are introduced since the drilling degree of freedom is used for the present shell element.^27,28

The objective of this research work is the development of a finite element for geometrically nonlinear shells with emphasis on the generality, robustness, and efficiency. The proposed element is a nine-node Lagrangian-type degenerated shell element. The element is based on the refined first-order shear deformation theory, and therefore is independent of the need of any shear correction factors. In the element formulation, an appropriate strain interpolation is employed to avoid both the shear and the membrane locking phenomena. Also, a suitable scheme is pursued for the expression of the rotational degrees of freedom so that the limitation on magnitude of the rotational increment in each load step can be removed. A full Gauss integration scheme is applied for the evaluation of the torsional stiffness.

The proposed formulation is computationally efficient and the test results showed good agreement. More specifically, the numerical examples of annular plate presented herein show good validity, efficiency, and accuracy for the developed nonlinear shell element.

Geometry and kinematics of the shell element

In Figure 1, the geometry of a nine-noded shell element with 6 degrees of freedom is illustrated.

Figure 1.

Geometry of nine-node shell element with 6 degrees of freedom.

The position vector and the displacement vector of the refined first-order shear deformation theory of nine-node shells can be simplified into ²⁴

$$\mathbf{Q}({\xi }_i)=\overline{\mathbf{Q}}({\xi }_j)+{\xi }_3(1-\frac{{\xi }_3^2}{3})\overline{\mathbf{V}}({\xi }_j)\mathbf{u}({\xi }_i)=\overline{\mathbf{u}}({\xi }_j)+{\xi }_3(1-\frac{{\xi }_3^2}{3})\overline{\mathbf{e}}({\xi }_j)$$

(1)

where $\mathbf{Q}$ and $\mathbf{u}$ denote the position vector and displacement vector of a generic point in the shell element, respectively; $\overline{\mathbf{Q}}$ and $\overline{\mathbf{V}}$ are the position vector of a point in the mid-surface and a normal vector to the mid-surface, respectively; $\overline{\mathbf{u}}$ and $\overline{\mathbf{e}}$ are the translational displacement vector and the fiber displacement vector of a point in the mid-surface, respectively.

Using the tensor transformation relationship of Green strain tensor and the natural strain, the transverse natural shear strains can obtain ²⁴

$${\tilde{E}}_{\alpha 3}^{\mathrm{so}}=\frac{1}{2}[\frac{\partial {\mathbf{Q}}_I}{\partial {\xi }_{\alpha }}{\overline{\mathbf{e}}}_I+\frac{\partial {\mathbf{u}}_J}{\partial {\xi }_{\alpha }}{\overline{\mathbf{V}}}_J+\frac{\partial {\mathbf{u}}_K}{\partial {\xi }_{\alpha }}{\overline{\mathbf{e}}}_K](1-{\xi }_3^2)$$

(2)

When the effective transverse shear energy ( $U_s$ ) is equal to the average transverse shear energy ( ${\overline{U}}_s$ ), the nominal uniform transverse shear strain ${\overline{E}}_{\alpha 3}$ is

$${\overline{E}}_{\alpha 3}=\frac{\int {\tilde{S}}_{\alpha 3}^s{\tilde{E}}_{\alpha 3}^{s}d{\xi }_3}{\int {\tilde{S}}_{\alpha 3}^{s}d{\xi }_3}=\frac{\int L_{\tau }({\xi }_3)L_{\gamma }({\xi }_3)d{\xi }_3}{\int L_{\tau }({\xi }_3)d{\xi }_3}{\tilde{E}}_{\alpha 3}^e$$

(3)

where $L_{\tau }({\xi }_3)$ and $L_{\gamma }({\xi }_3)$ are the distribution shape function of transverse shear stress and strain, respectively. ²⁶ Therefore, the ratio of transverse shear strain effective magnitude ${\tilde{E}}_{\alpha 3}^e$ to the nominal uniform shear strain ${\overline{E}}_{\alpha 3}$ is

$$f_{\xi }=\frac{{\tilde{E}}_{\alpha 3}^e}{{\overline{E}}_{\alpha 3}}=\frac{\int L_{\tau }({\xi }_3)d{\xi }_3}{\int L_{\tau }({\xi }_3)L_{\gamma }({\xi }_3)d{\xi }_3}$$

(4)

The refined first-order shear deformation theory can then be presented as a function of $(1-{\xi }_3^2)$ and $f_{\xi }$ factorizing the transverse shear strains in the traditional first-order shear deformation theory (equation (2)) as follows

$$\begin{array}{cc}\multicolumn{1}{c}{{\tilde{E}}_{\alpha 3}^{\mathrm{sn}}}& ={\tilde{E}}_{\alpha 3}^e·L_{\gamma }({\xi }_3)={\overline{E}}_{\alpha 3}·L_{\gamma }({\xi }_3)·f_{\xi }\hfill \\ \multicolumn{1}{c}{}& =\frac{1}{2}[\frac{\partial {\mathbf{Q}}_I}{\partial {\xi }_{\alpha }}{\overline{\mathbf{e}}}_I+\frac{\partial {\mathbf{u}}_J}{\partial {\xi }_{\alpha }}{\overline{\mathbf{V}}}_J+\frac{\partial {\mathbf{u}}_K}{\partial {\xi }_{\alpha }}{\overline{\mathbf{e}}}_K](1-{\xi }_3^2)f_{\xi }\hfill \end{array}$$

(5)

The limitation of modified first-order shear deformation theory ²⁴ is overcome using equation (5) for the correct analysis of thin and thick shells. Equation (5) is different from the usual natural strain relations of first-order shear deformation theory and the results in parabolic through-thickness distribution for the transverse shear strains.

The incremental displacement field within the element can be expressed as

$$\Delta \mathbf{u}({\xi }_i)={\displaystyle \sum _{a=1}^9}{N}_a({\xi }_j)\Delta {\mathbf{u}}^a({\xi }_3),i=1,2,3,j=1,2$$

(6)

where $N_a$ donates the two-dimensional quadratic Lagrangian interpolation function associated with node $a$ , $\Delta {\mathbf{u}}^a({\xi }_3)=[I_{3\times 3}|{\xi }_3{V}^a]\Delta U^a$ , and $\Delta U^a=\{\Delta {\overline{u}}_1^a,\Delta {\overline{u}}_2^a,\Delta {\overline{u}}_3^a,\Delta {\theta }_1^a,\Delta {\theta }_2^a,\Delta {\theta }_3^a\}$ . ²⁴

The shell element displays resultant membrane forces, moments, and transverse shear forces obtained by integration of stresses through the thickness using the equivalent constitutive equations

$$\begin{array}{cc}\multicolumn{1}{c}{}& \left\{\begin{array}{c}\hfill M_{\alpha \beta }^{(0)}\hfill \\ \multicolumn{1}{c}{M_{\alpha \beta }^{(1)}}\end{array}\right\}=\left[\begin{array}{cc}\hfill D_{\alpha \beta \gamma \delta }^{(0)}\hfill & \hfill D_{\alpha \beta \gamma \delta }^{(1)}\hfill \\ \multicolumn{1}{c}{D_{\alpha \beta \gamma \delta }^{(1)}}& \hfill D_{\alpha \beta \gamma \delta }^{(2)}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\tilde{E}}_{\gamma \delta }^m\hfill \\ \multicolumn{1}{c}{{\tilde{E}}_{\gamma \delta }^b}\end{array}\right\}\hfill \\ \multicolumn{1}{c}{}& \hfill \end{array}$$

(7a)

$$\left\{\begin{array}{c}\hfill Q_{\alpha 3}^{(0)}\hfill \\ \multicolumn{1}{c}{Q_{\alpha 3}^{(1)}}\end{array}\right\}=\left[\begin{array}{cc}\hfill D_{\alpha 3\beta 3}^{(0)}\hfill & \hfill D_{\alpha 3\beta 3}^{(2)}\hfill \\ \multicolumn{1}{c}{D_{\alpha 3\beta 3}^{(2)}}& \hfill D_{\alpha 3\beta 3}^{(4)}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\tilde{E}}_{\beta 3}^{\mathrm{sn}1}\hfill \\ \multicolumn{1}{c}{{\tilde{E}}_{\beta 3}^{\mathrm{sn}2}}\end{array}\right\}$$

(7b)

where $D_{\alpha \beta \gamma \delta }^{(i)}={\int }_{-h/2}^{h/2}{\tilde{C}}_{\alpha \beta \gamma \delta }(1,{\xi }_3,{\xi }_3^2)d{\xi }_3$ and $D_{\alpha 3\beta 3}^{(i)}={\int }_{-h/2}^{h/2}{\tilde{C}}_{\alpha 3\beta 3}(1,{\xi }_3^2,{\xi }_3^4)d{\xi }_3$ in which ${\tilde{C}}_{\alpha \beta \gamma \delta }$ is the fourth-order material tensor; $D_{\alpha \beta \gamma \delta }^{(i)}(i=0,1,2)$ are the extensional, bending–extensional coupling, and bending stiffness coefficients; and $D_{\alpha 3\gamma 3}^{(0)}$ and $D_{\alpha 3\gamma 3}^{(i)}(i=2,4)$ are the transverse shear and higher order transverse shear stiffness coefficients, respectively.

Assumed natural strains

Ma and Kanok-Nukulchai ²⁷ suggested the procedure for classifying the optimal locations of the sampling points to avoid locking phenomena. It was utilized by Han et al. ²⁴ The membrane and transverse shear strain fields are interpolated with the following sampling points in Figure 2.

Figure 2.

Sampling points for assumed strains of (a) membrane normal and transverse shear and (b) membrane shear.

Static equilibrium equation

The static equilibrium equation can be derived as

$${(f_{j\beta }{\tilde{S}}_{\alpha \beta })}_{,\alpha }+\tilde{\rho }{b}_j=0$$

(8)

where ${\tilde{S}}_{\alpha \beta }$ is the second and symmetric natural stress tensor, $b_j$ is the body force, and $f_{j\beta }$ is the deformation gradient which is defined as

$$f_{j\beta }=\frac{\partial x_j}{\partial {\xi }_{\beta }}$$

(9)

and $\tilde{\rho }$ is the mass density of the body defined as

$$\tilde{\rho }={\rho }_0{\tilde{J}}_t/J$$

(10)

in which ${\tilde{J}}_t$ is $|\partial x_i/\partial {\xi }_{\alpha }|$ , ${\rho }_0$ is the mass density of the undeformed body, and $J$ is $|\partial P_I/\partial x_j|$ .

The weak form of the static problem can be constructed by Galerkin’s weighted residual method as

$$\begin{array}{cc}\multicolumn{1}{c}{G^e(\mathbf{u}\mathbf{,}\delta \mathbf{u})=}& -\underset{B^e}{\int }[{({\tilde{S}}_{\alpha \beta }{f}_{l\beta })}_{,\alpha }+\tilde{\rho }{b}_l]\delta u_l\mathrm{dV}\hfill \\ \multicolumn{1}{c}{}& +\underset{\partial B^e}{\int }[{\tilde{S}}_{\alpha \beta }{f}_{l\beta }{n}_{\alpha }-T_l^{(N)}]\delta u_l\mathrm{dA}=0\hfill \end{array}$$

(11)

where $\mathbf{u}$ and $\delta \mathbf{u}$ are the displacement field and weight field, respectively; $B^e$ and $\delta B^e$ are the volume and its boundary surface of the parental element, respectively. $T_l^{(N)}$ is a natural traction vector associated with area element. Applying the Gauss–Green theorem, the above equation can be rewritten in the canonical form as

$$\underset{B^e}{\int }{\tilde{S}}_{\alpha \beta }{f}_{l\beta }\delta u_{l,\alpha }\mathrm{dV}-\underset{B^e}{\int }\tilde{\rho }{b}_l\delta u_l\mathrm{dV}-\underset{\partial B^e}{\int }{T}_l^{(N)}\delta u_l\mathrm{dA}=0$$

(12)

The weight field of Galerkin’s method has the same function space as the displacement field, so the weight field can be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{\delta \mathbf{u}({\xi }_1,{\xi }_2,{\xi }_3)}& ={\displaystyle \sum _{b=1}^9}\frac{\partial \mathbf{x}({\xi }_1,{\xi }_2,{\xi }_3)}{\partial {\mathbf{U}}^b}\delta {\mathbf{U}}^b\hfill \\ \multicolumn{1}{c}{}& ={\displaystyle \sum _{b=1}^9}{N}_b({\xi }_j)\left[{\mathbf{I}}_{3\times 3}|{\xi }_3{\mathbf{V}}^b\right]\delta {\mathbf{U}}^b\hfill \end{array}$$

(13)

where $\delta {\mathbf{U}}^b$ denotes the virtual nodal displacements. Substituting equation (13) into equation (12) yields

$${\displaystyle \sum _{b=1}^9}\delta U_j^b(K_j^b-R_j^b)=0$$

(14)

in which $K_j^b$ is the internal force vector and $R_j^b$ is the generalized force vector.

Since $\delta U_j^b$ can be arbitrary, equation (14) reduces to

$$K_j^i-R_j^i=0\mathrm{or}\mathbf{K}\mathbf{(}\mathbf{U}\mathbf{)}-\mathbf{R}=0$$

(15)

The drilling degree of freedom will be linked to the in-plane twisting mode of the mid-surface by a penalty functional. ²⁹

Tangent stiffness and arc-length control method

To solve the nonlinear equation by the Newton–Raphson method, the linearization is required, and the result is the element tangent stiffness matrix in the form

$$\partial {\mathbf{K}\mathbf{(}\mathbf{U}}_m^i\mathbf{)}\Delta {\mathbf{U}}_m^i={\mathbf{R}}_m-{\mathbf{K}\mathbf{(}\mathbf{U}}_m^i\mathbf{)}$$

(16)

where $\partial \mathbf{K}$ is the tangent stiffness with the value of a displacement ${\mathbf{U}}_m^i$ in which the subscript $m$ refers to the current mth load step and superscript $i$ refers to the current ith iteration step. The linearization of internal force vector is required, and the result is the tangent stiffness in the form

$$\partial {\mathbf{K}}_{\mathrm{ji}}^{\mathrm{ba}}=\frac{\partial K_j^b}{\partial U_i^a}=\frac{\partial }{\partial U_i^a}\left(\underset{B^e}{\int }{\tilde{S}}_{\alpha \beta }\frac{\partial x_l}{\partial {\xi }_{\beta }}\frac{\partial }{\partial {\xi }_{\alpha }}\left(\frac{\partial x_l}{\partial U_j^b}\right)\mathrm{dV}\right)$$

(17)

In order to trace the full path of load–displacement, it is necessary to use the arc-length control method, since the postbuckling study involves a highly distorted load–displacement path. With consideration of the incremental load step, equation (15) can be rewritten in the following form which has (n + 1) unknowns

$$G\mathbf{(}\mathbf{U},\lambda \mathbf{)}=\mathbf{K}\mathbf{(}\mathbf{U}\mathbf{)}-\lambda \mathbf{F}=0$$

(18)

where $\mathbf{K}\mathbf{(}\mathbf{U}\mathbf{)}$ is the internal force vector which is generally a nonlinear function of $\mathbf{U}$ , $\lambda$ is a loading intensity, and $\mathbf{F}$ is a reference load vector. The nonlinear equations are augmented by an additional constraint equation (19) ^30,31

$$\begin{array}{cc}\multicolumn{1}{c}{f\mathbf{(}\mathbf{U},\lambda \mathbf{)}:}& {\displaystyle \sum _{k=1}^n}{\chi }_k{(u_k^m-u_k^{m-1})}^2\hfill \\ \multicolumn{1}{c}{}& +{\chi }_{n+1}{\alpha }^2{({\lambda }^m-{\lambda }^{m-1})}^2\hfill \\ \multicolumn{1}{c}{}& ={\displaystyle \sum _{k=1}^n}{\chi }_k\Delta u_k^2+{\chi }_{n+1}{\alpha }^2\Delta {\lambda }^2=c^2\hfill \end{array}$$

(19)

where $\alpha$ denotes a load coefficient and $c$ is the prescribed arc-length. For arc-length control method, ${\displaystyle \sum _{k=1}^n}\Delta u_k^2+{\alpha }^2\Delta {\lambda }^2=c^2$ when all ${\chi }_k=1$ .

Linearization of equation (18) and the constraint equation $f\mathbf{(}\mathbf{U},\lambda \mathbf{)}$ with respect to $\mathbf{U}$ and $\lambda$ about the previous solution ${\mathbf{(}\mathbf{U}}^{i,m},{\lambda }^{i,m}\mathbf{)}$ leads to the linearized equilibrium and constraint equation. Finally, the nonlinear equation augmented by a constraint equation can be written in matrix form

$$\left[\begin{array}{cc}\hfill \partial {\mathbf{K}}^{i,m}\hfill & \hfill -\mathbf{F}\hfill \\ \multicolumn{1}{c}{C_1^{\mathrm{T}}}& \hfill C_2\hfill \end{array}\right]\left\{\begin{array}{c}\hfill \Delta {\mathbf{U}}^{i,m}\hfill \\ \multicolumn{1}{c}{\Delta {\lambda }^{i,m}}\end{array}\right\}=\left\{\begin{array}{c}\hfill {\mathbf{R}}^{i,m}\hfill \\ \multicolumn{1}{c}{r^{i,m}}\end{array}\right\}$$

(20)

where $\partial {\mathbf{K}}^{i,m}$ is the tangent stiffness for a known displacement vector ${\mathbf{U}}^{i,m}$ , $C_1^{\mathrm{T}}=\{2{\chi }_1(u_1^{i,m}-u_1^{m-1}),\dots ,2{\chi }_n(u_n^{i,m}-u_n^{m-1})\}$ , $C_2=2{\chi }_{n+1}{\alpha }^2({\lambda }^{i,m}-{\lambda }^{m-1})$ , ${\mathbf{R}}^{i,m}$ is the unbalanced force vector, and superscript i,m denotes the quantity at the ith iteration of the mth solution step. The incremental quantities ${\mathbf{U}}^{i+1,m}={\mathbf{U}}^{i,m}+\Delta {\mathbf{U}}^{i,m}$ and ${\lambda }^{i+1,m}={\lambda }^{i,m}+\Delta {\lambda }^{i,m}$ .

Numerical examples

Several numerical examples are solved to test the performance of the shell element in geometrically nonlinear applications. Examples included isotropic plates to check some crucial features and shells for the comparisons and further developments. Before the various sections, nonlinear examples of hemispherical shells are treated. First, nonlinear isotropic hemispherical shell with $18\circ$ open and full hemispherical shell examples are presented for comparison with previously published reference. Second, isotropic twisted beam example is presented for geometrically linear and nonlinear comparisons. Finally, more general shell examples treated with various types for parameter studies.

Nonlinear analysis of hemispherical shell

The problem configuration of hemispherical shell with an $18\circ$ circular cutout at its pole is defined in Figure 3. This example is well known in the literature. The nonlinear version of this problem should be more difficult to fully demonstrate that the formulation can undergo large displacements. Here, one quarter of the shell is modeled using $8\times 8$ , $12\times 12$ , and $16\times 16$ nine-node shell elements. The loads are increased by a factor of 425 to compare with the results by Sze et al. ³² The good performance of the large displacement analysis is shown in Figure 4. The average number of iterations per load step is 6 for a total of 20 load steps. The initial and deformed configurations of the hemispherical shell with an 18° circular cutout at its pole are depicted in Figure 5.

Figure 3.

Definition of loading points and initial configuration of ${18}^{\circ }$ open hemispherical shell.

Figure 4.

Load–deflection curves for $18\circ$ open hemispherical shell.

Figure 5.

Initial configuration and deformed configuration of $18\circ$ open hemispherical shell: (a) initial configuration and (b) deformed configuration.

The next example considered is the full pinched hemispherical shell under two inward and two outward opposite forces (Figure 6). In this case, we consider a full hemisphere with no hole. The radial displacement at points A and B is monitored in Figure 7. Comparisons of the present results with those of Arciniega and Reddy ³³ are in close agreement, but the reference results are slightly stiffer. Arciniega and Reddy ³³ used a regular mesh of $2\times 2Q81$ elements. In this study, the mesh is composed of three parts, with either $4\times 4$ or $8\times 8$ elements in each (Figure 6). The average number of iterations per load step is 6.7 for a total of 31 load steps. The initial and deformed configurations of the full pinched hemisphere are presented in Figure 8.

Figure 6.

Definition of loading points and initial configuration of hemispherical shell.

Figure 7.

Load–deflection curves for a full hemispherical shell.

Figure 8.

Initial configuration and deformed configuration of full hemispherical shell: (a) initial configuration and (b) deformed configuration.

The third example considered is the hemispherical shell with an 18° circular cutout at its pole. In the first example, the thickness $h=0.04$ is used. The four times smaller thickness, $h=0.01$ , is used in the nonlinear tests. For the larger radius-to-thickness ratio $R/h=1000$ , the shell elements are more prone to locking. The displacement in the direction of the force at the point where this force is applied is reported. We see that there are no large differences between the Wisniewski and Turska ³⁴ and the present shell element, but the reference results are slightly stiffer. Wisniewski and Turska ³⁴ used a regular mesh of $16\times 16$ HW47 elements. The solution curves for the inward displacement under the force are shown in Figure 9.

Figure 9.

Load–deflection curves for $18\circ$ open hemispherical shell ( $h=0.01$ ) at point B (inward displacement).

Nonlinear analysis of twisted beam

The initial geometry of the beam is twisted, but the initial strain is equal to 0. The beam is clamped at one end and loaded by a unit force at the other. The force is applied either along the y-axis (in-plane) or along the x-axis (out-of-plane), see Figure 10. Here, we provide the results for the latter case only. The data are as follows: $E=2.9\times {10}^7$ , $\nu =0.22$ , the length $L=12$ , the width $W=1.1$ , and the twist is $90\circ$ . This example belongs to the set of tests of MacNeal and Harder, ³⁵ and flat shell elements have problems with this example. In our computations, the $3\times 18$ element mesh and a very small thickness $h=0.0032$ are used. The results of linear analysis for $F_0=1\times {10}^{-6}$ are given in Table 1, where the displacement in the direction of force at end point is presented. Although the shell is extremely thin, the results are not corrupted by the membrane locking and are quite close to the reference value. The nonlinear load–deflection curves obtained are shown in Figures 11 and 12 where the displacement in the direction at end point is monitored.

Figure 10.

Definition of loading points and initial configuration of twisted beam.

Table 1.

Displacement in the direction of force at end point.

Element	HW47 Wisniewski and Turska 34	EADG5A Wisniewski and Turska 34	Belytschko et al. 36	Present (mesh)
Displacement	1.2895	1.2882	1.2940	1.2927 (3 × 18)
				1.2918 (2 × 12)

Figure 11.

Load–deflection curves for a twisted beam under in-plane load.

Figure 12.

Load–deflection curves for a twisted beam under out-of-plane load.

Nonlinear analysis of C-shape narrow plate

The C-shape narrow plate problems are shown in Figure 13 with sequence of deformed configuration with tip bending moment. The total length of C-shape narrow plate is $L=120+120+120$ , width is $b=10$ , Young’s modulus is $E=1000$ , Poisson’s ratio is $v=0.0$ , and the thickness is $t=1.0$ . As shown in Figure 13, the tip bending moment $M_A$ is loaded at the free end, and deformed shapes are illustrated in Figure 13. The load parameter of tip bending moment ranges from 0 to $6M_{A}L/(7\pi \mathrm{EI})$ .

Figure 13.

Sequence of deformed finite element meshes corresponding to the increasing loads of C-shape narrow plate. End B is fixed. End A is free with bending moment.

Nonlinear analysis of pinched cylinder of a clamped cylinder

Another well-known benchmark problem of finite deformation is the pinched cylindrical shell under a point load (Figure 14(a)). This problem has been investigated by Balah and Al-Ghamedy, ³⁹ Sze et al., ³² and Arciniega and Reddy. ³³ The material properties are as follows: Young’s modulus $E=2.0685\times {10}^7$ and Poisson’s ratio $v=0.3$ . The cylinder length is 3.048 and the radius is 1.016 with thickness of 0.03. One-quarter of the structure is analyzed with $14\times 14$ and $16\times 16$ element meshes. The load–deflection curve is shown in Figure 15 and compared with Sze et al. ³² In Figure 15, it can be observed that the load–deflection curve goes beyond the highest physically possible deflection of the loading point. The initial and deformed configurations of the pinched cylinder are depicted in Figure 14(b).

Figure 14.

Problem definition and a deformed configuration of pinched cylinder with free edge: (a) initial configuration and (b) deformed configuration.

Figure 15.

Load–deflection curve of clamped cylindrical shell under pinching point load.

Nonlinear analysis of slit annular plate

Annular plate under end out-of-plane shear force

In order to validate the present shell element, 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 (Figure 16). This beautiful benchmark problem was considered by Sansour and Kollmann ³⁷ and Arciniega and Reddy, ³³ among others, and Sze et al. ³² presented ABAQUS results for the problem. In particular, we use the ABAQUS results presented in a recent article of Sze et al. ³² for popular benchmark problems of nonlinear shell analysis (because of the tabulated data). The geometry and elastic material properties for the isotropic case are given by

$$E=2.1\times {10}^7,\nu =0.0,R_i=6.0,R_o=10.0,h=0.03$$

(21)

Figure 16.

Slit annular plate loaded with various end forces F.

for a maximum distributed force of $F=1.45$ . A regular mesh of $5\times 40$ elements is considered in this analysis. Solutions obtained by the mesh of $2\times 16$ elements are too stiff and suffer from locking. Arciniega and Reddy ³³ used 80 load steps with equal load steps of 0.01 to the $\mathrm{F}=0.8$ , but the computation is performed by the Newton–Raphson method with 54 load steps (with equal load steps of 0.015) in this study. The standard convergence tolerance for the residual forces is set to be $1\times {10}^{-6}$ . The shear load versus displacement curves for two characteristic points are depicted in Figures 17 and 18. Solutions obtained with the present formulation are in good agreement with those obtained by Sze et al. ³² Figure 19 shows the deformed configuration of the isotropic annular plate for load level $\mathrm{F}=1.45$ . It is clear that the plate undergoes large displacements at the corresponding loading of F = 1.45.

Figure 17.

Load–deflection curves for the slit annular plate lifted by line force F.

Figure 18.

Load–deflection curves for the slit annular plate lifted by line force F.

Figure 19.

Deformed configuration of the slit annular plate (5 × 40 mesh) for F = 1.45.

Annular plate under end moment

We examine an annular clamped plate subjected to end loads (Figure 20). The load versus displacement curves for two characteristic points are depicted in Figure 21. Figure 22 shows the deformed configuration of the isotropic annular plate for load level $\mathrm{F}=10.0$ .

Figure 20.

Slit annular plate loaded with various end forces F.

Figure 21.

Load–deflection curves for the slit annular plate deformed by two point forces.

Figure 22.

Deformed configuration of the slit annular plate (5 × 40 mesh) for F = 10.0.

Annular plate under end in-plane force ( $F_y(+)$ )

We examine an annular clamped plate subjected to end loads (Figure 23). To activate the postbuckling equilibrium path we have to include a small perturbation in the system. Otherwise, only in-plane displacements in the plate will occur (which is called the fundamental equilibrium path beyond which the critical load will be unstable). For the secondary path, we assume that the annular plate is presented with an imperfect displacement of $R_0/{\delta }_z=100$ at loading edge. That is, in addition to the in-plane compressive load there is a small out of plane loading that causes the postbuckling response of the annular plate. The tip displacement of the annular plate versus the axial load is shown in Figure 24. It is clearly observed that a well-defined pre-buckling regime ( $F<F_{\mathrm{cr}}$ ) and postbuckling regime ( $F>F_{\mathrm{cr}}$ ) exist. Figure 25 shows the deformed configuration of the isotropic annular plate for load level $\mathrm{F}=0.38$ .

Figure 23.

Slit annular plate loaded with various end forces F.

Figure 24.

Load–deflection curves for the slit annular plate lifted by line force F.

Figure 25.

Deformed configuration of the slit annular plate (5 × 40 mesh) for F = 0.38.

Annular plate under end in-plane force ( $F_y(-)$ )

We examine an annular clamped plate subjected to end loads (Figure 26). The load versus displacement curves for two characteristic points are depicted in Figure 27. Figure 28 shows the deformed configuration of the isotropic annular plate for load level $\mathrm{F}=0.55$ .

Figure 26.

Slit annular plate loaded with various end forces F.

Figure 27.

Load–deflection curves for the slit annular plate lifted by line force F.

Figure 28.

Deformed configuration of the slit annular plate (5 × 40 mesh) for F = 0.55.

Annular plate under end in-plane shear force ( $F_x(+)$ )

We examine an annular clamped plate subjected to end loads (Figure 29). The load versus displacement curves for two characteristic points are depicted in Figure 30. Figure 31 shows the deformed configuration of the isotropic annular plate for load level $\mathrm{F}=1.40$ .

Figure 29.

Slit annular plate loaded with various end forces F.

Figure 30.

Load–deflection curves for the slit annular plate lifted by line force F.

Figure 31.

Deformed configuration of the slit annular plate (5 × 40 mesh) for F = 1.40.

Annular plate under end in-plane shear force ( $F_x(-)$ )

We examine an annular clamped plate subjected to end loads (Figure 32). The load versus displacement curves for two characteristic points are depicted in Figure 33. Figure 34 shows the deformed configuration of the isotropic annular plate for load level $\mathrm{F}=1.40$ .

Figure 32.

Slit annular plate loaded with various end forces F.

Figure 33.

Load–deflection curves for the slit annular plate lifted by line force F.

Figure 34.

Deformed configuration of the slit annular plate (5 × 40 mesh) for F = 0.2875.

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 $3\times 3$ or $6\times 6$ meshes, as shown in Figure 35. Three different thickness values are considered, that is, h = 3.175, 6.35, and 12.7 mm. The problem has been considered in Sabir and Lock, ⁴⁰ Horrigmoe and Bergan, ⁴¹ Rhiu and Lee, ⁴² Sze et al., ³² and Arciniega and Reddy. ³³ The arc-length control method is used to trace the equilibrium path, and a tolerance coefficient of $1.0\times {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=3102.75{\mathrm{N}/\mathrm{mm}}^2$ , and Poisson’s ratio $\nu =0.3$ .

Figure 35.

Hinged cylindrical shell; curved edges are free, straight edges are hinged (P = 1000).

For the case of h = 12.7 mm, the average number of iterations is 5, and the curve of the central deflection versus load is given in Figure 36, showing the good correlation between the present solution and the reference solutions.

Figure 36.

Load–deflection curves for the hinged shell (h = 12.7 mm).

For the case of h = 6.35 mm, the solution derived with $24\times 24$ four-node shell elements by Sze et al. ³² is used as the major reference. The average number of iterations in each load step is 4. It is seen that the solution in this study is very close to the solution of Sze et al. ³² in Figure 37.

Figure 37.

Load–deflection curves for the hinged shell (h = 6.35 mm).

In the third case of this analysis, even a further reduction in the panel thickness was presumed taking h = 3.175 mm, what resulted in a rather complicated shape of the equilibrium path. The curve of the central deflection versus load is given in Figure 38.

Figure 38.

Load–deflection curves for the hinged shell (h = 3.175 mm).

Nonlinear analysis of pullout of an open-ended cylindrical shell

The example consists of a cylindrical shell with open ends, which is pulled by a pair of radial forces. The geometrical data and material properties are presented in Figure 39. The problem has been considered in Sansour and Kollmann, ³⁷ Sze et al., ³² and Arciniega and Reddy, ³³ among others. Owing to symmetry, one-eighth of the shell is modeled and a mesh with six elements along the circumferential direction and four elements along the longitudinal direction is adopted. Points A, B, C, and D are monitored in terms of radial displacement in Figure 40. The displacements for points A–D are compared with the results obtained by Sze et al. ³² It is noticed that the response of the shell has 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 deformed configuration under the load level of P = 40,000 is portrayed in Figure 41.

Figure 39.

Open-end cylindrical shell subjected to radial pulling forces (P = 40,000).

Figure 40.

Load–deflection curves for the open-end cylindrical shell.

Figure 41.

Deformed configuration of the cylindrical shell (6 × 4 mesh) for P = 40,000.

Nonlinear analysis of pinched cylindrical shell

The pinched cylindrical shell mounted on rigid end diaphragms is shown in Figure 42. The example consists of a cylindrical shell subjected to two radial opposite forces. The dimensions and the material properties for one-eighth of the shell are illustrated in Figure 42. Due to symmetry, one-eighth of the shell is modeled and the employed mesh for nine-node shell elements is $16\times 16,20\times 20,24\times 24$ , and $30\times 30$ . This test was carried out in Sansour and Kollmann, ³⁷ El-Abbasi and Meguid, ⁴³ and Sze et al., ³² among others. Figure 43 shows the load–displacement curves. The result of the present shell element gives a very close agreement with that of the reference as shown in Sze et al. ³² It is noted that the finite element discretization can slightly change the nature of the response for coarse meshes so that local snap-through appears. In order to compute the minimum mesh refinement that captures the nature of the fine-mesh solution, further study is required. Adaptive meshing technique may be one way to alleviate the difficulties of computation of the minimum mesh refinement. In Figure 44, we show the deformed configuration of the isotropic cylindrical shell for load level ${\mathrm{P}}_1=0.997\mathrm{P}$ .

Figure 42.

Pinched cylindrical shell mounted on rigid end diaphragms (P = 12,000).

Figure 43.

Load–deflection curves for the pinched cylindrical shell.

Figure 44.

Deformed configuration of the pinched cylindrical shell (24 × 24 mesh) for ${\mathrm{P}}_1=0.997\mathrm{P}$ .

Conclusion

The finite element formulation of an improved nine-node shell element is presented for the solution of geometrically nonlinear problems of annular plates, cylindrical, and spherical shells. Geometrically nonlinear analysis of shells using an ELF is proposed. To overcome the limit in the very thin plates and shells, shear and membrane locking is avoided using assumed natural strains. The stress values are required in the formation of the internal force vector and the element tangent stiffness matrix. To avoid the adverse effect of inaccurate stresses on the finite element solution, better stress evaluation method should be adopted. To achieve this goal, strain interpolation scheme similar to the one used in the computation of linear stiffness matrix is first employed for the strain evaluation. Then, based on the interpolated strain field, stresses at integration points are calculated. The numerical results available show that due to such a treatment, the element performance in solving nonlinear problems with distorted mesh can be greatly improved. The present shell element formulation is based on refined first-order transverse shear deformation theory, which results in parabolic through-thickness distribution of the transverse shear strains and stresses. It eliminates the need for shear correction factors in the first-order theory. The scheme for finite rotation of shell normal is used to remove the limitation on magnitude of the rotational increment in each load step in geometrically nonlinear analysis. A full Gauss integration scheme is applied for the evaluation of the membrane, bending, shear, and torsional stiffness. The present solutions show very good agreement with the other numerical solutions. Especially, the annular plate and cylindrical shells may provide the benchmark to test the results of geometrically nonlinear analysis for shell structures.

It should be pointed out that the development of an “absolutely” good element is impossible and one has to be satisfied with those possessing certain desirable properties. It can be expected that the study reported will contribute to further investigation of different approaches for finite element improvement, and thereby to the development of more effective finite elements for engineering applications.