The discrete model and the analysis of a spherical shell by finite equilibrium elements

S. Kalanta*, J. Atkočiūnas**, T. Ulitinas*** *Vilnius Gediminas Technical University, Saulėtekio av. 11, 10223 Vilnius, Lithuania, E-mail: stanislovas.kalanta@vgtu.lt **Vilnius Gediminas Technical University, Saulėtekio av. 11, 10223 Vilnius, Lithuania, E-mail: juozas.atkociunas@vgtu.lt ***Vilnius Gediminas Technical University, Saulėtekio av. 11, 10223 Vilnius, Lithuania, E-mail: tomas.ulitinas@vgtu.lt


Introduction
Computer-aided analysis of structures requires the development of their discrete model with the finite number of the degrees of freedom.Usually, the finite element method is used because it is universal, effective, and easy for evaluating the boundary conditions.There are three major modifications of the finite element method [1], including the methods, when only the function of displacement or the function of the internal forces (equilibrium elements) is approximated, and a composite method.In the equilibrium finite element method, equations of statics are satisfied completely or almost completely, while the geometric equations of interrelations are satisfied approximately [2,3].The problem of calculating the internal forces, displacements and strains under the given load is reduced to deriving the equations of the displacement method.These equations are derived by summing up the stiffness values of the elements according to the algorithm, regulated by equations of statics.
In this paper, the discretization of a spherical shell is thoroughly investigated by the equilibrium circular finite elements, which enables the authors to apply the unified methodology to the analysis of the elastic and elasticplastic shells [2][3][4][5][6].The equilibrium finite element is symmetrically loaded for flat spherical shells.The bending moments and axial forces are described by the second and first-degree polynomials.The element's differential statics equations, describing the balance between the internal and external forces, are replaced with algebraic equilibrium equations presented in the Bubnov-Galerkin method.The mathematical model and the calculation algorithm of the internal forces and displacements in the shell analysis are developed and formulated, using the equations of statics and geometry.The analysis and description of the method of shell discretization by equilibrium finite elements are still lacking in the scientific literature.

General data on the discrete shell model
The schematic views of the spherical shell and the internal forces are given in Figs. 1 and 2. The stress state of the shell under the symmetrical load is defined by the vector of the internal forces while the external distributed load vector is expressed as follows The bending moments M ρ (ρ), M φ (ρ) and the axial forces of the intensity functions N ρ (ρ), N φ (ρ), are the stress vector's components (their positive directions are shown in Fig. 2).
Fig. 1 The spherical shell Fig. 2 The internal forces of the spherical shell The vectors-functions S(ρ) and p(ρ) are related by the differential equations of statics [ ] ( ) ( ) where the differential operator is expressed as A where R 0 is the radius of curvature of the shell.Each element of the circular shell is marked by the index k = 1, 2,…, r (r is the number of the elements).The nodal (calculation) points in the element are indexed by i = 1, 2, 3 as shown in Fig. 3.The shell is considered in the cylindrical (ρ, φ, z) coordinate system, while the element is analysed in the local coordinate system ξ (Fig. 4).
All the elements of the shell are connected by the boundary nodes in the main nodes 3, 5, 7 of the discrete model (Fig. 3).The second-order circular element is used for discretization, while the internal forces are approximated by the second-degree polynomials, which are created in [2].The discrete shell model is regular for circular elements of the same width.The load can be distributed over the surface of the finite elements or concentrated in the main node.
Fig. 3 Discretization of the flat spherical shell by four finite circular elements The assumption is made that physical properties of the material (the elastic modulus E and Poisson's ratio ν), shell thickness t and the intensity of the distributed load in the element are constant.
The functions of internal forces, defining the stress state, are approximated by all the vector-functions of the internal forces having a finite number of elements ( ) ( ) where [H k (ξ)] is the interpolation matrix of the element's internal forces, developed in the local coordinate system ξ is presented in the second and fifth rows of Table 1.S k is the vector of the internal forces of the element's nodes.The interpolation points of the internal forces are the nodal points of the finite element.The unknown coefficients of the function are the components of the vector S k .Thus, the stress state of the discrete model is defined by the internal forces' vector S = (S 1 , S 2 ,…, S k ,…, S r ) T .It is one of the unknowns of the computational shell problem.A generalized nodal force vector P k and its dual displacement vector u k of the nodal point are constructed for the element k.The number of components of each vector m k defines the element's degree of freedom.The work of the element's internal forces S k must be equal to the work of the external forces P k .The local forces (distributed in the unit of area) or the concentrated forces acting in the main element's nodes may be considered to be the generalized forces.In the first case, the vector u k is composed of the integral displacements of the whole element, while, in the second case, these are the local node displacements.Thus, the global displacement vector u = (u 1 , u 2 ,…, u m ) T describes the deformation state of the discrete model of the shell.Here, m is the degree of freedom of the discrete model.The relationship between the global displacements u and the local displacements u k and is expressed by the equation where [B k ] is displacement compatibility matrix of the k-th element of the shell.Generalized forces and displacements should be selected so that the equations of statics for the directions of the displacements u are valid for the internal and boundary nodes of all elements.
Given the basic data on the equilibrium finite elements, we can obtain the main element's dependences and develop a mathematical model for calculating the circular shell's element.

The dependences and matrices of the finite circular element
Bending moments of the shell are approximated by the second-degree polynomials, and the axial force -by the linear functions.The internal forces at the nodal points of the finite element are shown in Fig. 4. The vector S k of the internal forces is presented in the first row of Table 1.The interpolation matrix of the internal forces [H k (ξ k )] is given elsewhere in the Table 1.The components ( ) , , ; where ρ k2 is the coordinate of the second node in the global coordinate system (ρ, φ, z); 2b k is the width of the finite element.Algebraic internal equilibrium equations of the finite element are obtained by inserting the interpolation functions of the internal forces Eq. ( 1) into the Eq. ( 2) and differentiating ( ) The operator of the algebraic equations [A k (ξ k )] is presented in Table 2.The operator [A k (ξ k )] depends on the coordinates ξ k .Equations of statics of the element are derived, using Bubnov-Galerkin method, which states that the equilibrium is guaranteed at some nodal points of linear independent statics equations.In the considered case, the element's nodes are considered to be the collocation nodes.The same influence functions are taken for both equations and the matrix of the influence function is constructed Equations of statics of the element Eq. ( 5) are obtained by using the formula The matrix [A ek ] of the equations of statics of the element is given in the rows 4-7 in Table 3.It is argued that the intensity of the load distributed over the element's surface is constant, while the surface area of the element is equal to the horizontal projection area because a flat shell is considered.

Table 2
The algebraic operator of the element's equilibrium equations The vector of the external node forces is equivalent to the distributed element's load , where the coefficient of the first node (i = 1) is c = 1, while the coefficient of the third node (i = 3) is c = -1.

Sub-matrices [A k1 ], [A ek ] and [A k3
] make the matrix of equations of statics of the shell element The matrices [A k1 ] and [A k3 ] define the relationship between the internal forces of the element's boundary nodes and the generalized forces P k1 and P k3 .They are given in the first and the last three rows in Table 3.The expressions of the forces P k3 and P k10 are obtained by using the dependence Their coefficients are given in the rows 3 and 10 in Table 3.
The geometric equations of the element are as follows 0 where the flexibility matrix is obtained by the formula The matrix of the infinitely small flexibility element is is the Poisson's ratio; t k is the element thickness.Flexibility matrix of the shell element is presented in Table 4.
The matrices [A k ] and [D k ] are used for developing the stiffness matrices of the spherical shell elements.It is clear that they depend on the position of the elements and, therefore, are constructed individually for each element.For this purpose, only the values of the element's width b k , the second node coordinate  k2 and physical parameters should be inserted into the obtained expressions.

The mathematical model of the problem of shell analysis
The mathematical model of the problem of the analysis of the elastic shell's displacements and internal forces consists of m algebraic equations of statics and n geometric equations.The unknowns are the n-dimensional vector S of internal forces and the m-th dimensional displacement vector u.The system of Eqs. ( 9), (10) defines the stress and strain state of the construction.The main system of equations of equilibrium finite elements is a connectinglink to the elastic-plastic shells calculation [7][8][9][10].The development of the matrix [A] of the equilibrium equations' coefficients and flexibility matrix [D] is briefly discussed below.Eqs. ( 5) of statics of the elements k = 1, 2,…, r and the main shell's nodes of the equilibrium equations make the discrete model of statics Eq. ( 9) of the shell.
Equations of statics of the discrete model's node j, where the elements k and l meet (Fig. 5), consist of the equilibrium equations of the bending moments and axial and shear forces where F n,j is the intensity of the normally distributed load in the circular element  j .Equations of statics of the shear forces are derived, using the dependence (7).Geometric Eq. ( 10) are formed by geometric Eq. ( 8) of all finite elements, formulated, taking into account the compatibility equations of displacements Eq. (3).

The flexibility matrix
obtained from the mathematical model Eqs.( 9), ( 10) by eliminating the internal forces where [K] is the global stiffness matrix of the construction.
Displacements are calculated by the formula uF .The internal forces in the elements' nodes are calculated by the formula (12) or Though stiffness matrix [K] can be calculated by the formulas (13), in the developed program, it is constructed in a usual way, based on the finite element stiffness matrices and using the algorithm described in the article [3].

The analysis of the element accuracy and convergence
A firmly fixed shell, subjected to normally distributed loading of the intensity p = 1 kN/m 2 , and having the radius of curvature R 0 = 1 m is considered (Fig. 6).The values of the bending moments of the multiplex mesh do not differ from the values obtained with 4 finite elements.The values of the axial forces (up to the factor 0 pR ) are presented in Table 5.We can see that, when eight finite elements are taken, sufficiently accurate calculation results are obtained.Two axial forces N φ at the middle shell point 5 (Fig. 7), corresponding to the nodes of two adjacent ele-ments' nodes are shown in column 4 of Table 5.We can see that, with the increase of the number of elements, the discontinuities are getting smaller.The discontinuities are quite small, when 8 elements are taken.
Computational analysis allows the authors to conclude that the created shell element is sufficiently accurate.It also confirms the statement about the higher accuracy of the equilibrium finite elements compared to the displacement elements [5].Therefore, the created element can be effectively used in the elastic-plastic shell analysis and optimization [11][12][13].

Conclusion
1.The presented element dependencies allow the equations, describing nodal displacements of a discrete model, to be directly derived by using the stiffness matrix of the elements (similar to the method of the displacement elements).They are formulated according to the algorithm described in the paper [3] and using the flexibility matrix of the element presented in Table 4.
2. The created element can be effectively used for the elastic-plastic spherical shell analysis as well as for formulating and solving the optimization problems.
3. The performed computational analysis, using the mesh of the elements of various density, has shown that the accuracy and convergence of the calculation results are high.This is particularly important for the analysis of the elastic-plastic shells and for solving the optimizationnonlinear programming problems, whose solution success largely depends on their size (the number of elements).

THE DISCRETE MODEL AND THE ANALYSIS OF A SPHERICAL SHELL BY FINITE EQUILIBRIUM ELEMENTS S u m m a r y
The paper presents the equilibrium finite element discretization of symmetrically loaded spherical flat shells.It is based on Castigliano principle.A new second-order equilibrium finite element is suggested, and the equilibrium and physical equations, obtained for it by using the Bubnov-Galiorkin method, are presented.A mathematical model for solving the problem of the elastic shell computation is created, based on the above equations.The methodology is illustrated by a numerical example.The results are obtained, using a computer-aided program developed by the authors.The calculation results, obtained using the mesh of the elements of various density, show that the accuracy of the created element and the convergence of the results are high.

Fig. 4
Fig. 4 The layout of the finite element of the shell The relationship between the global coordinate ρ k and the local coordinate ξ k is defined by the dependencies is the modulus of elasticity of the element k; k 

Fig. 5
Fig. 5 Forces acting on the main node of the shell Statics equations matrix     T A B A   

Fig. 6 A
Fig. 6 A computational scheme of the spherical shell The shell is divided into four circular elements of the same width.A discrete model fragment is shown in Fig. 7.The values of nodal displacements are given up to the factor E pR / 0

Fig. 7
Fig. 7 Displacement and internal forces of shell (up to the factors E pR / 0

Table 1 Interpolation
matrix of internal forces ( )

Table 3
Statics equation matrix of the spherical shell element

Table 5
The values of the axial forces