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.

Keywords

Geometrically nonlinear analysis element-based Lagrangian shell element modified first-order shear deformation theory assumed natural strain

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.

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

Q (ξ_{i}) = \bar{Q} (ξ_{j}) + ξ_{3} (1 - \frac{ξ_{3}^{2}}{3}) \bar{V} (ξ_{j}) u (ξ_{i}) = \bar{u} (ξ_{j}) + ξ_{3} (1 - \frac{ξ_{3}^{2}}{3}) \bar{e} (ξ_{j})

(1)

where $Q$ and $u$ denote the position vector and displacement vector of a generic point in the shell element, respectively; $\bar{Q}$ and $\bar{V}$ are the position vector of a point in the mid-surface and a normal vector to the mid-surface, respectively; $\bar{u}$ and $\bar{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}}_{α 3}^{so} = \frac{1}{2} [\frac{\partial Q_{I}}{\partial ξ_{α}} {\bar{e}}_{I} + \frac{\partial u_{J}}{\partial ξ_{α}} {\bar{V}}_{J} + \frac{\partial u_{K}}{\partial ξ_{α}} {\bar{e}}_{K}] (1 - ξ_{3}^{2})

(2)

When the effective transverse shear energy ( $U_{s}$ ) is equal to the average transverse shear energy ( ${\bar{U}}_{s}$ ), the nominal uniform transverse shear strain ${\bar{E}}_{α 3}$ is

{\bar{E}}_{α 3} = \frac{\int {\tilde{S}}_{α 3}^{s} {\tilde{E}}_{α 3}^{s} d ξ_{3}}{\int {\tilde{S}}_{α 3}^{s} d ξ_{3}} = \frac{\int L_{τ} (ξ_{3}) L_{γ} (ξ_{3}) d ξ_{3}}{\int L_{τ} (ξ_{3}) d ξ_{3}} {\tilde{E}}_{α 3}^{e}

(3)

where $L_{τ} (ξ_{3})$ and $L_{γ} (ξ_{3})$ are the distribution shape function of transverse shear stress and strain, respectively.²⁶ Therefore, the ratio of transverse shear strain effective magnitude ${\tilde{E}}_{α 3}^{e}$ to the nominal uniform shear strain ${\bar{E}}_{α 3}$ is

f_{ξ} = \frac{{\tilde{E}}_{α 3}^{e}}{{\bar{E}}_{α 3}} = \frac{\int L_{τ} (ξ_{3}) d ξ_{3}}{\int L_{τ} (ξ_{3}) L_{γ} (ξ_{3}) d ξ_{3}}

(4)

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

\begin{matrix} {\tilde{E}}_{α 3}^{sn} & = {\tilde{E}}_{α 3}^{e} \cdot L_{γ} (ξ_{3}) = {\bar{E}}_{α 3} \cdot L_{γ} (ξ_{3}) \cdot f_{ξ} \\ = \frac{1}{2} [\frac{\partial Q_{I}}{\partial ξ_{α}} {\bar{e}}_{I} + \frac{\partial u_{J}}{\partial ξ_{α}} {\bar{V}}_{J} + \frac{\partial u_{K}}{\partial ξ_{α}} {\bar{e}}_{K}] (1 - ξ_{3}^{2}) f_{ξ} \end{matrix}

(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

Δ u (ξ_{i}) = \sum_{a = 1}^{9} N_{a} (ξ_{j}) Δ u^{a} (ξ_{3}), i = 1, 2, 3, j = 1, 2

(6)

where $N_{a}$ donates the two-dimensional quadratic Lagrangian interpolation function associated with node $a$ , $Δ u^{a} (ξ_{3}) = [I_{3 \times 3} | ξ_{3} V^{a}] Δ U^{a}$ , and $Δ U^{a} = {Δ {\bar{u}}_{1}^{a}, Δ {\bar{u}}_{2}^{a}, Δ {\bar{u}}_{3}^{a}, Δ θ_{1}^{a}, Δ θ_{2}^{a}, Δ θ_{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{matrix} {\begin{matrix} M_{α β}^{(0)} \\ M_{α β}^{(1)} \end{matrix}} = [\begin{matrix} D_{α β γ δ}^{(0)} & D_{α β γ δ}^{(1)} \\ D_{α β γ δ}^{(1)} & D_{α β γ δ}^{(2)} \end{matrix}] {\begin{matrix} {\tilde{E}}_{γ δ}^{m} \\ {\tilde{E}}_{γ δ}^{b} \end{matrix}} \end{matrix}

(7a)

{\begin{matrix} Q_{α 3}^{(0)} \\ Q_{α 3}^{(1)} \end{matrix}} = [\begin{matrix} D_{α 3 β 3}^{(0)} & D_{α 3 β 3}^{(2)} \\ D_{α 3 β 3}^{(2)} & D_{α 3 β 3}^{(4)} \end{matrix}] {\begin{matrix} {\tilde{E}}_{β 3}^{sn 1} \\ {\tilde{E}}_{β 3}^{sn 2} \end{matrix}}

(7b)

where $D_{α β γ δ}^{(i)} = \int_{- h / 2}^{h / 2} {\tilde{C}}_{α β γ δ} (1, ξ_{3}, ξ_{3}^{2}) d ξ_{3}$ and $D_{α 3 β 3}^{(i)} = \int_{- h / 2}^{h / 2} {\tilde{C}}_{α 3 β 3} (1, ξ_{3}^{2}, ξ_{3}^{4}) d ξ_{3}$ in which ${\tilde{C}}_{α β γ δ}$ is the fourth-order material tensor; $D_{α β γ δ}^{(i)} (i = 0, 1, 2)$ are the extensional, bending–extensional coupling, and bending stiffness coefficients; and $D_{α 3 γ 3}^{(0)}$ and $D_{α 3 γ 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.

Static equilibrium equation

The static equilibrium equation can be derived as

{(f_{j β} {\tilde{S}}_{α β})}_{, α} + \tilde{ρ} b_{j} = 0

(8)

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

f_{j β} = \frac{\partial x_{j}}{\partial ξ_{β}}

(9)

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

\tilde{ρ} = ρ_{0} {\tilde{J}}_{t} / J

(10)

in which ${\tilde{J}}_{t}$ is $| \partial x_{i} / \partial ξ_{α} |$ , $ρ_{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{matrix} G^{e} (u, δ u) = & - \int_{B^{e}} [{({\tilde{S}}_{α β} f_{l β})}_{, α} + \tilde{ρ} b_{l}] δ u_{l} dV \\ + \int_{\partial B^{e}} [{\tilde{S}}_{α β} f_{l β} n_{α} - T_{l}^{(N)}] δ u_{l} dA = 0 \end{matrix}

(11)

where $u$ and $δ u$ are the displacement field and weight field, respectively; $B^{e}$ and $δ 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

\int_{B^{e}} {\tilde{S}}_{α β} f_{l β} δ u_{l, α} dV - \int_{B^{e}} \tilde{ρ} b_{l} δ u_{l} dV - \int_{\partial B^{e}} T_{l}^{(N)} δ u_{l} 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{matrix} δ u (ξ_{1}, ξ_{2}, ξ_{3}) & = \sum_{b = 1}^{9} \frac{\partial x (ξ_{1}, ξ_{2}, ξ_{3})}{\partial U^{b}} δ U^{b} \\ = \sum_{b = 1}^{9} N_{b} (ξ_{j}) [I_{3 \times 3} | ξ_{3} V^{b}] δ U^{b} \end{matrix}

(13)

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

\sum_{b = 1}^{9} δ 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 $δ U_{j}^{b}$ can be arbitrary, equation (14) reduces to

K_{j}^{i} - R_{j}^{i} = 0 or K (U) - 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 {K (U}_{m}^{i}) Δ U_{m}^{i} = R_{m} - {K (U}_{m}^{i})

(16)

where $\partial K$ is the tangent stiffness with the value of a displacement $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 K_{ji}^{ba} = \frac{\partial K_{j}^{b}}{\partial U_{i}^{a}} = \frac{\partial}{\partial U_{i}^{a}} (\int_{B^{e}} {\tilde{S}}_{α β} \frac{\partial x_{l}}{\partial ξ_{β}} \frac{\partial}{\partial ξ_{α}} (\frac{\partial x_{l}}{\partial U_{j}^{b}}) dV)

(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 (U, λ) = K (U) - λ F = 0

(18)

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

\begin{matrix} f (U, λ) : & \sum_{k = 1}^{n} χ_{k} {(u_{k}^{m} - u_{k}^{m - 1})}^{2} \\ + χ_{n + 1} α^{2} {(λ^{m} - λ^{m - 1})}^{2} \\ = \sum_{k = 1}^{n} χ_{k} Δ u_{k}^{2} + χ_{n + 1} α^{2} Δ λ^{2} = c^{2} \end{matrix}

(19)

where $α$ denotes a load coefficient and $c$ is the prescribed arc-length. For arc-length control method, $\sum_{k = 1}^{n} Δ u_{k}^{2} + α^{2} Δ λ^{2} = c^{2}$ when all $χ_{k} = 1$ .

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

[\begin{matrix} \partial K^{i, m} & - F \\ C_{1}^{T} & C_{2} \end{matrix}] {\begin{matrix} Δ U^{i, m} \\ Δ λ^{i, m} \end{matrix}} = {\begin{matrix} R^{i, m} \\ r^{i, m} \end{matrix}}

(20)

where $\partial K^{i, m}$ is the tangent stiffness for a known displacement vector $U^{i, m}$ , $C_{1}^{T} = {2 χ_{1} (u_{1}^{i, m} - u_{1}^{m - 1}), \dots, 2 χ_{n} (u_{n}^{i, m} - u_{n}^{m - 1})}$ , $C_{2} = 2 χ_{n + 1} α^{2} (λ^{i, m} - λ^{m - 1})$ , $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 $U^{i + 1, m} = U^{i, m} + Δ U^{i, m}$ and $λ^{i + 1, m} = λ^{i, m} + Δ λ^{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.

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 2 Q 81$ 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.

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.

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}$ , $ν = 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.

Table 1.

Displacement in the direction of force at end point.

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

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 $6 M_{A} L / (7 π EI)$ .

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

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}, ν = 0.0, R_{i} = 6.0, R_{o} = 10.0, h = 0.03

(21)

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 $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 $F = 1.45$ . It is clear that the plate undergoes large displacements at the corresponding loading of 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 $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} / δ_{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_{cr}$ ) and postbuckling regime ( $F > F_{cr}$ ) exist. Figure 25 shows the deformed configuration of the isotropic annular plate for load level $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 $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 $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 $F = 1.40$ .

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 {N / mm}^{2}$ , and Poisson’s ratio $ν = 0.3$ .

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.

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.

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.

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.

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 $P_{1} = 0.997 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.

Footnotes

Acknowledgements

The authors would like to express profound gratitude to Prof. W. Kanok-Nukulchai in AIT for his encouragement while preparing this article.

Academic Editor: Stephane PA Bordas

Declaration of conflicting interests

The authors declare that there is no conflict of interests regarding the publication of this article.

Funding

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (MEST) (No. 2011-0028531).

References

Ahmad

Irons

Zienkiewicz

. Analysis of thick and thin shell structures by curved finite elements. Int J Numer Meth Eng 1970; 2: 419–451.

Huang

Hinton

. A new nine node degenerated shell element with enhanced membrane and shear interpolation. Int J Numer Meth Eng 1986; 22: 73–92.

MacNeal

. Derivation of element stiffness matrices by assumed strain distributions. Nucl Eng Des 1982; 33: 1049–1058.

Dvorkin

Bathe

. A continuum mechanics based four-node shell element for general non-linear analysis. Eng Computation 1984; 1: 77–88.

Kim

Lomboy

Han

. A co-rotational 8-node assumed strain shell element for postbuckling analysis of laminated composite plates and shells. Comput Mech 2003; 30: 330–342.

Bathe

Dvorkin

. A formulation of general shell elements—the use of mixed interpolation of tensorial components. Int J Numer Meth Eng 1986; 22: 697–722.

Nguyen-Xuan

Rabczuk

Bordas

. A smoothed finite element method for plate analysis. Comput Meth Appl Mech Eng 2008; 197: 1184–1203.

Nguyen-Thanh

Rabczuk

Nguyen-Xuan

. A smoothed finite element method for shell analysis. Comput Meth Appl Mech Eng 2008; 198: 165–177.

Rodrigues

Natarajan

Ferreira

AJM

. Analysis of composite plates through cell-based smoothed finite element and 4-noded mixed interpolation of tensorial components techniques. Comput Struct 2014; 135: 83–87.

10.

Kanok-Nukulchai

Wong

. Element-based Lagrangian formulation for large deformation nonlinear analysis. Comput Struct 1988; 30: 967–974.

11.

Zhang

Kim

. Linear and geometrically nonlinear analysis of plates and shells by a new refined non-conforming triangular plate/shell element. Comput Mech 2005; 36: 331–342.

12.

Wisniewski

. Finite rotation shells. Basic equations and finite elements for Reissner kinematics. Barcelona: CIMNE Springer, 2010.

13.

Wisniewski

Turska

. Improved four-node Hellinger-Reissner elements based on skew coordinates. Int J Numer Meth Eng 2008; 76: 798–836.

14.

Wisniewski

Turska

. Improved four-node Hu-Washizu elements based on skew coordinates. Comput Struct 2009; 87: 407–424.

15.

Firl

Bletzinger

. Shape optimization of thin walled structures governed by geometrically nonlinear mechanics. Comput Meth Appl Mech Eng 2012; 237–240: 107–117.

16.

Shin

Lee

. Development of a strain-smoothed three-node triangular flat shell element with drilling degrees of freedom. Finite Elem Anal Des 2014; 86: 71–80.

17.

Kim

Liu

Han

. A resultant 8-node solid-shell element for geometrically nonlinear analysis. Comput Mech 2005; 35: 315–331.

18.

Ahmed

van der Meer

Sluys

. A geometrically nonlinear discontinuous solid-like shell element (DSLS) for thin shell structures. Comput Meth Appl Mech Eng 2012; 201–204: 191–207.

19.

Pandya

Kant

. Higher-order shear deformable theories for flexure of sandwich plates finite element evaluations. Int J Solids Struct 1988; 24: 1267–1286.

20.

Reddy

. A general nonlinear 3rd-order theory of plates with moderate thickness. Int J Nonlin Mech 1990; 25: 677–686.

21.

Noor

Burton

Bert

. Computational models for sandwich panels and shells. Appl Mech Rev 1996; 49: 155–199.

22.

Rolfes

Rohwer

. Improved transverse shear stress in composite finite element based on first order shear deformation theory. Int J Numer Meth Eng 1997; 40: 51–60.

23.

Tanov

Tabiei

. A simple correction to the first-order shear deformation shell finite element formulations. Finite Elem Anal Des 2000; 35: 189–197.

24.

Han

Tabiei

Park

. Geometrically nonlinear analysis of laminated composite thin shells using a modified first-order shear deformable element-based Lagrangian shell element. Compos Struct 2008; 82: 465–474.

25.

Han

Kanok-Nukulchai

Park

. A refined finite element for first-order plate and shell analysis. Struct Eng Mech 2011; 40: 191–213.

26.

Jung

Han

. Shear buckling responses of laminated composite shells using a modified 8-node ANS shell element. Compos Struct 2014; 109: 119–129.

27.

Kanok-Nukulchai

. On the application of assumed strained methods. In: Kanok-Nukulchai

Karasudhi

(eds) Structural engineering and construction, achievement, trends and challenges. Bangkok, Thailand: AIT, 1989, pp.1168–1175.

28.

Han

Kim

Kanok-Nukulchai

. An element-based 9-node resultant shell element for large deformation analysis of laminated composite plates and shells. Struct Eng Mech 2004; 18: 807–829.

29.

Kanok-Nukulchai

. A simple and efficient finite element for general shell analysis. Int J Numer Meth Eng 1979; 14: 179–200.

30.

Lee

Kanok-Nukulchai

. A nine-node assumed strain finite element for large deformation analysis of laminated shells. Int J Numer Meth Eng 1998; 42: 777–798.

31.

Chaiomphob

Kanok-Nukulchai

Nishino

. An automatic arc-length control algorithm for tracing equilibrium paths of nonlinear structures. Struct Eng Earthquake Eng JSCE 1988; 5: 205–208.

32.

Sze

Liu

. Popular benchmark problems for geometric nonlinear analysis of shells. Finite Elem Anal Des 2004; 40: 1551–1569.

33.

Arciniega

Reddy

. Tensor-based finite element formulation for geometrically nonlinear analysis of shell structures. Comput Meth Appl Mech Eng 2007; 196: 1048–1073.

34.

Wisniewski

Turska

. Four-node mixed Hu–Washizu shell element with drilling rotation. Int J Numer Meth Eng 2012; 90: 506–536.

35.

MacNeal

Harder

. A proposed standard set of problems to test finite element accuracy. Finite Elem Anal Des 1985; 1: 3–20.

36.

Belytschko

Stolarski

Liu

. Stress projection for membrane and shear locking in shell finite elements. Comput Meth Appl Mech Eng 1985; 51: 221–258.

37.

Sansour

Kollmann

. Families of 4-node and 9-node finite elements for a finite deformation shell theory. An assessment of hybrid stress, hybrid strain and enhanced strain elements. Comput Mech 2000; 24: 435–447.

38.

Panasz

Wisniewski

. Remarks on the selective reduced integration in 9-node shell elements. In: Pietraszkiewicz

Kreja

(eds) Shell structures: theory and applications. London: Taylor & Francis, 2010, pp. 273–276.

39.

Balah

Al-Ghamedy

. Finite element formulation of a third-order laminated finite rotation shell element. Comput Struct 2002; 80: 1975–1990.

40.

Sabir

Lock

. The application of finite elements to the large deflection geometrically non-linear behaviour of cylindrical shells. In: Brebia

Tottenham

(eds) Variational methods in engineering. Southampton: Southampton University Press, 1973, pp.7/66–7/75.

41.

Horrigmoe

Bergan

. Nonlinear analysis of free form shells by flat finite elements. Comput Meth Appl Mech Eng 1978; 16: 11–35.

42.

Rhiu

Lee

. A nine-node finite element for analysis of geometrically non-linear shells. Int J Numer Meth Eng 1988; 24: 581–604.

43.

El-Abbasi

Meguid

. A new shell element accounting for through-thickness deformation. Comput Meth Appl Mech Eng 2000; 189: 841–862.

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

Abstract

Keywords

Introduction

Geometry and kinematics of the shell element

Assumed natural strains

Static equilibrium equation

Tangent stiffness and arc-length control method

Numerical examples

Nonlinear analysis of hemispherical shell

Nonlinear analysis of twisted beam

Nonlinear analysis of C-shape narrow plate

Nonlinear analysis of pinched cylinder of a clamped cylinder

Nonlinear analysis of slit annular plate

Annular plate under end out-of-plane shear force

Annular plate under end moment

Annular plate under end in-plane force ( F y ( + ) )

Annular plate under end in-plane force ( F y ( − ) )

Annular plate under end in-plane shear force ( F x ( + ) )

Annular plate under end in-plane shear force ( F x ( − ) )

Nonlinear analysis of hinged shell

Nonlinear analysis of pullout of an open-ended cylindrical shell

Nonlinear analysis of pinched cylindrical shell

Conclusion

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

References

Annular plate under end in-plane force ( $F_{y} (+)$ )

Annular plate under end in-plane force ( $F_{y} (-)$ )

Annular plate under end in-plane shear force ( $F_{x} (+)$ )

Annular plate under end in-plane shear force ( $F_{x} (-)$ )