Nonlinear stability analysis of shallow cylindrical shells with initial imperfections

The present article describes numerical procedures for geometrically nonlinear stability analysis of shallow flexible cylindrical shells. Computing models are based on the theory of plates with account of transverse shear deformations. The finite difference energy method of discretization is used for reducing the continuum problem to finite dimensional problem. Solution procedures for discrete nonlinear problem are based on initial value formulation. The solution of Cauchy problem is performed using the self-correcting Runge-Kutta method with the constrained equation. For stability analysis of shells with initial imperfections small asymmetric perturbations are introduced in the initial conditions of the Cauchy problem for finding asymmetric shapes of shell equilibrium.


Introduction
In the general case the solution of the stability problem is reduced to the identification of qualitative changes in a behavior of the system at changes in the structure of the system.The system is structurally stable if sufficiently small perturbations lead to correspondingly small changes in its behavior.A measure of the deviation from the unperturbed states is internal state parameters (generalized coordinates).
The stability analysis of a shell structure studies its behavior under small perturbation of parameters including in the energy functional, as well as the initial shape, boundaries and boundary conditions.Such a generalization of the concept of stability is important for practical calculations as in real structures it is possible not only imperfections in the form or a small unbalance in the load but also a heterogeneity of the physical properties of material and a deviation from set boundaries and boundary conditions.In the analysis of shells the definition of such initial imperfections in the form that lead to the largest drop of the critical load compared to the same load for the perfect system becomes important.It allows to obtain the lower estimation of the actual critical load of a real shell.relationships between the strains and the displacements have the form [1][2][3] where e 11 , e 22 , e 12linear and angular deformations; κ 11 , κ 22 , κ 12the changes of curvatures and curvature of torsion; e 13 , e 23the transverse shear deformations; u, v, wdisplacements in the direction of coordinate axes x, y and z; θ 1 , θ 2angles of rotation of shell cross-sections.
The relationships (1) constructed with consideration of transverse shear deformations which allows to analyze the shells with medium thickness and made from the materials with low shear stiffness.If we assume that e 13 =e 23 =0, from the expressions (1) the geometric relations of Kirchhoff-Love theory of shells and plates are followed [4,5].
For analyze the stress-strain state and stability of cylindrical shells the finite difference energy method is used [6,7].Lagrange functional of the theory of shells has the form: where εstrains vector; Dthe elasticity matrix linking stresses and strains; qthe external load vector whose components correspond to the five components of the displacement vector; uthe vector of nodal displacements; Ωthe area occupied by the shell.
In this method a grid is applied on the region Ω of the changing of independent variables.The desired function u, delivers the stationary value to the functional (2), approximately determined by its values at the nodes, and the derivatives of u are replaced by finite differences.Using finite-difference formulas for displacements and its derivatives and replace integration by summation over the area Ω occupied by the plate, the functional (2) reduces to its discrete analog, which is a scalar function of a vector argument.Euler equations for the Lagrange functional (2) are the equations of equilibrium in displacements.
The equilibrium equations of the shell in the finite difference energy method can be represented as: the gradient of the potential energy of deformation; Qthe normalized vector of nodal loads; pthe parameter of the load.
All the components of the displacement vector u and the load parameter p is assumed to dependent on some continuation parameter s.Each value s will uniquely correspond to certain equilibrium s p s , u . A continuous sequence of equilibrium states of the system is the equilibrium curve.For its determination method of differentiating with respect to parameter is used [8,9].It brings the system of nonlinear equations (3) to Cauchy problem for systems of ordinary differential equations with the constrained equation: The solution of Cauchy problem is performed using the self-correcting Runge-Kutta method [1,10]: where a positive scalar parameter.The components of the gradient vector and Hessian matrix elements are calculated using the formulas [1,10] based on expressions (1) and ( 2): The application of procedure of the self-correcting Runge-Kutta method for solution of the system (3) provides an opportunity to build the efficient numerical algorithm with high precision and the absence of accumulation of the calculations errors.Self-correcting scheme of the Runge-Kutta method on the step m includes [1,10] The solution of equations ( 8) is performed by Newton-Raphson method.The calculations are performed for one iteration.The expression (7), (8)  the solution obtained using the method of differentiating with respect to parameter is specified at a certain step using Newton-Raphson procedure.

Results
The developed method for studying the stability of the equilibrium forms of the shell structures were used in the analysis of shallow cylindrical shells under transverse uniformly distributed load with intensity q z (Fig. 1).Boundary conditions are hinged support along the contour.The load is uniformly distributed load with intensity E p h a q q z / / 2 k is made by the method of differentiating with respect to parameter (4).The arc length of the equilibrium curve is taken as a parameter of continuation.The numerical solution was obtained using a self-correcting Runge -Kutta method (5-9) with increment 5 , 0 's .In Fig. 2 equilibrium curves of the considered panels with symmetric and asymmetric deformation are shown (p is load parameter, W c is the normal displacement in the center of the shall).Curve 1 illustrates the symmetric deformation of the shell.

Fig. 1. The general form of a shallow cylindrical shell
When solving this problem using the Crisfield scheme (4) it is necessary to separate limit points from the bifurcation points, and thus to determine the sign of the increment of continuation parameter.On each step of leading parameter the value of stiffness S p =U T Q is calculated.Vector U is determined from the solution of the equation . The moment of the simultaneous change of the signs of the matrix determinant )) ( ( 2m k s W u and the parameter S p indicates the transition trough the limit point.In the analysis of shallow shells the definition of such initial imperfections that cause the largest drop of the critical load compared to the same load for the perfect system becomes important.It allows us to obtain the lower estimation of the actual critical load of a real shell.Numerical studies have shown that the most dangerous are the initial imperfections, which are determined by its eigenform, obtained from the solution of the eigenvalue problem.
Small asymmetric perturbations corresponding to the first asymmetric eigenforms are where βsmall number, about 10 -2 h; i, jthe number of half waves in the direction of the x and y axes respectively.Given small initial asymmetric imperfections (10) in the initial stages does not affect the stress-strain state practically but when approaching to the bifurcation point it leads to the transition to the asymmetric curve.In this case, the critical points corresponding to the asymmetric shape of the stability loss, defined not as the bifurcation points but as limit points.
For the transition to the asymmetric branch of the equilibrium curve 2 (Fig. 2) it is given the values 1 i , 2 j in the formula (10), which corresponds a small asymmetry with respect to the line x=a.For the transition to the asymmetric branch of the equilibrium curve 3 (Fig. 2) it is given the values 1 i , 1 j in the formula (10), which corresponds a small asymmetry with respect to the two straight lines x=a and y=a.
The correction of the first asymmetric branch of equilibrium curve in the vicinity of the bifurcation point is performed as the following: after entering on the asymmetric branch the initial imperfection is deleted and the reverse motion along the equilibrium curve by setting the negative increment parameter is implemented.This method let specify the limit load in asymmetric deformation of shells.

Conclusions
The analysis of the obtained results allows to draw some conclusions about the applied computational schemes: 1. Rapid restructuring of the equilibrium forms in some levels of loading makes a significant error in the solution of the problem using the standard procedures of Runge-Kutta method.Self-correcting form of Runge-Kutta method (5-9) compensates the error in the solution from the previous step which allows to determine the entire equilibrium curve with a sufficient degree of accuracy.2. Numerical integration of Cauchy problem (4) with the constraint equation, describing the spherical surface, requires a much smaller integration step than applying the constrained equation, describing the plane.If the calculations are performed with a large step, then in the second case there are oscillations of the obtained solutions with respect to the true equilibrium position.
3. Comparison of techniques based on the self-correcting method of increments and Runge-Kutta method shows that to obtain results with the accuracy within 1% is required to set the step for self-correcting method of increments four times less than for Runge-Kutta method.It leads to the same cost of computer time for the numerical problem solution of both methods.
4. Euler method has a slow convergence.Using of this method in numerical computations of nonlinear deformable shells is associated with high computer time due to the necessity of choosing a small step continuation parameter, especially in the stability problems.
5. The most effective methods for solving problems of this class are Cholesky decomposition and the LDL T -factorization.Cholesky decomposition, sometimes called the square root method is possible only for symmetric, positive definite matrix of the algebraic equations system.In the study of stability of thin-walled structures in the supercritical region Hessian matrix of potential energy of deformation becomes negative definite, therefore the use of Cholesky decomposition becomes impossible.When using LDL Tfactorization the condition of positive definiteness of the matrix is not imposed, so this procedure can be successfully applied for construction of the equilibrium curves.
S-P Seminar 2017, Theoretical Foundation of Civil Engineering 1 introduced in the initial conditions of Cauchy problem (4) for finding asymmetric shapes of shell equilibrium:

Fig. 2 .
Fig. 2. Equilibrium curves of a cylindrical shell with curvatures k 1 =60, k 2 =0 when symmetric and asymmetric deformation are the fourth order accurate, so step s ' can be chosen significantly larger than for Euler method which has the first order accuracy.Optionally,