Geometrically Nonlinear Behavior of Axisymmetric Thin Spherical Shells

This article presents a geometrically nonlinear formulation of a two node axisymmetric shell element. The geometrically nonlinear formulation is based on the Total Lagrangian approach and the material behavior is assumed to be linearly elastic. The spherical arc-length procedure is used to obtain the pre-buckling, buckling and post-buckling deformation path. Some numerical examples are solved to demonstrate geometrically nonlinear behavior of elastic thin spherical shells subjected to various loads.


Introduction
The analysis of shells always presents a challenge because of difficulties in formulation as well as numerical calculation. Moreover, shells of revolution have received considerable attention for applications such as underwater pressure hulls, space vehicles, pressure vessels, and pressure tanks. The most frequently investigated structures are cylindrical, conical, and spherical shells and their combinations [1,2]. The prediction of buckling loads and deformation paths of structural members is important in the design of various engineering components [3]. Many researches have been made on buckling of axisymmetric thin spherical shells under various loading and boundary conditions, such as, [4][5][6][7][8].
In the linear analysis, displacements and strains developed in a structure are small. Namely, the structure geometry remains unchanged during the loading process and therefore, linear strain approximations can be used. However, the structure geometry changes continuously during the loading process, and this fact is taken into account in geometrically nonlinear analysis. Generally, two types of buckling analyses methods which are linear and nonlinear are used in the finite element methods. Linear buckling analysis is also called eigenvalue buckling analysis because it gives rise to an algebraic eigenvalue problem. Eigenvalue buckling analysis method has two versions, classical and fully linearized buckling analysis. Linear buckling analysis gives good prediction of buckling loads if the prebuckling rotations are negligible. However, shell structures have considerable prebuckling rotations and linear or eigenvalue buckling analysis alone is not sufficient to predict the stability limit of these structures [9]. In order to correctly predict the load carrying capacity of a structure, the buckling and post buckling behavior of structures should be obtained.
Three Lagrangian kinematic descriptions are usually used for finite element analysis of geometrically nonlinear structures: The Total Lagrangian (TL), the Updated Lagrangian (UL) and the Co-Rotational (CR) formulations. Both the Total Lagrangian and the Updated Lagrangian formulations include all kinematic nonlinear effects due to large displacements, large rotations and large strains, but whether the large strain behavior is modeled approximately depends on the constitutive relations specified. The only advantage of using one formulation rather than the other lies in its greater numerical efficiency. The Co-Rotational formulation is the most recent of the three and the least developed one. Unlike the others, its domain of application is limited by a priori kinematic assumptions: displacements and rotations may be arbitrarily large, but deformations must be small. Because of this restriction, the CR has not penetrated the major general-purpose FEM codes that provide to nonlinear analysis [10][11][12].
In this study a geometrically nonlinear formulation of axisymmetric shells is given. The element has two nodes and three degrees of freedom per each node. The Total Lagrangian method, all static and kinematic variables at the current state are referred to the initial configuration, is used. To maintain slope and displacement continuity, the displacements within the element have to be uniquely determined by the nodal displacements and the position.
The incremental equilibrium equations are solved using the spherical arc-length method with automatic load increments. The choice of the incremental load factor sign in the predictor phase of arc-length methods is known to be of paramount importance in determining the success of such procedures in tracing unstable equilibrium paths [13]. Feng et al. (1996) have proposed a direction prediction criterion whereby the sign of the predictor load factor is made consistent with the sign of the internal product between the previous converged incremental displacement and the current tangential solution. This criterion is insensitive to the presence of bifurcations [14,15].
Some numerical examples are performed by a computer program which is written by the author in MATLAB software. The buckling and post buckling behaviors of elastic spherical shells subjected to different loads are demonstrated.

Element Geometry, Stresses and Strains
An axisymmetric shell element subjected to bending and membrane forces is given in Figure 1. The strain vector ε can be written using the Kirchhoff-Love assumptions which excludes transverse shear deformations, i. e. the angle ϕ does not vary, as where 0 ε and L ε are the linear and nonlinear strain vectors, respectively [16]. In this equation first two terms denote membrane strains and second two terms denote bending strains, respectively. For such elements the element geometry is described in terms of the coordinates of middle surface nodes.
The stress vector can be related to the strain vector using an elasticity matrix D For an isotropic shell the elasticity matrix is where E is elasticity modulus, ν is Poisson ratio and t is shell thickness [17].  Figure 2 shows a two-noded straight element which can be used to analyze the axisymmetric shells. The element has the radial and axial displacements, u and w , and a rotation, / dw ds β = , at each node. Then, the local displacements of a node i can be defined as The two-noded element has 6 dofs and they are determined by the element node displacements as To maintain slope and displacement continuity, the displacements within the element have to be uniquely determined by the nodal displacements a and the position s. Thus, the element displacements in local coordinates can be written as where N is the matrix that contains the shape functions. At the node i, the following equation can be written between local and global displacements sin cos (7) where Τ is the transformation matrix [18].

Stiffness Matrix
Equilibrium conditions between internal and external forces must be satisfied. Thus, the equilibrium equation of the nonlinear system can be written as where 0 B and L B represents linear and nonlinear straindisplacement matrices, respectively.
Taking the variation of Eq. (8) and using Eqs. (2 and 9), we have in which the matrix σ K is known as initial stress matrix.
And 0 K represents small displacement stiffness matrix, i. e., The matrix L K is known as large displacement matrix and is given by Consequently, the tangent stiffness matrix T K can be given as

Solution Method
The equilibrium equation of the nonlinear system can be written as where R , F and P represent the out-of-balance force vector, the internal force vector and the externally applied load vector respectively, and ( ) i λ ∆ is the load-level parameter In order to compute the nodal displacements, it is necessary to solve the system of nonlinear equilibrium equations using an incremental/iterative method. The load controlled Newton-Raphson method is the earliest method that is used to trace the equilibrium path. This method is based on the linearization of the equilibrium equations at a prescribed load level, that is, (14) is kept constant during iterations. The iterations are performed until the residual is smaller than a prescribed tolerance. This method can trace the load-displacement curve before the occurrence of a limit point, but generally it will fail to converge beyond this point.
In order to compute the nodal displacements, the spherical arc-length algorithm is used. The convergence criterion is chosen as -5 G G δu < 10 ∆u (15) where δ and ∆ parameters indicate iterative and incremental quantities, respectively. The choice of the sign of the incremental load factor in the predictor phase of arclength methods is significant for tracing unstable equilibrium paths. If the wrong sign is predicted, the solution sequence ''doubles back" on the original load-detection curve and the arc-length method fails to trace the complete path. Moreover, in the case of sharp snapback situations, the predictor criterion also becomes very important. It is seen that using the predictor criterions of Feng et al. [14] can overcome these problems successfully. Figure 3 shows the description of spherical cap which is subjected to the ring load with variable radius. The cap is clamped at its border. Load location is determined by load location ratio, e / θ θ = . In this ratio, θ and θ represent angles of the cap support and the ring loads according to the vertical line intersecting the apex, respectively. In this manner, load location ratio may take values from 0 to 1. As the ratio increases the ring load approaches to the cap support. Geometric and material properties are also given in the figure. The behaviors of the spherical cap for varying location of the ring load are given in Figures 4-9. In these figures, v A and v L denotes vertical deflections of the apex and the loaded node, respectively. To mesh the spherical cap 20 and 40 elements are used for e=0.00-0.80 and e=0.90, respectively. The element matrices and the load vectors are formed using 2x2 Gauss quadrature rule. It is clearly seen that the snap-through behavior begins when e=0.30 approximately and, in addition to snap-through behavior also snapback behavior occurs when location ratio close to e=0.50. Moreover, number of limit points increase as the load location ratios increase, after e=0.60 especially.

A Hinged Axisymmetric Shell Under Ring LOAD
The previous axisymmetric shell under ring load is analyzed for a different boundary condition, that is, the cap is assumed hinged at its border. The load-displacement curves of spherical cap for different location of the ring loads are given in Figures 10-15. Unlike the clamped condition, both snap-through behavior and snapback behavior occurs even e=0.00. Moreover, number of limit points increase and curves become very complex as the load location ratios increase. In these figures, also v A and v L denotes vertical deflections of the apex and the loaded node, respectively. To mesh the spherical cap 20 and 40 elements are used for e=0.00-0.80 and e=0.90, respectively.

Conclusion
In this paper, a geometrically nonlinear formulation based on the Total Lagrangian approach for a conical frustum element was given assuming the material behavior to be linearly elastic. The Arc-length method was used to analyze the post-buckling behavior. The behaviors of the spherical cap for varying location of the ring load and different boundary conditions were studied. Numerical examples showed that the snap-through and snapback behaviors of the cap vary depending on the load location. The complexity of the structure's behavior increases as the ring load moves from the apex to the cap support and also boundary condition of the cap effects the behavior significantly.