AN ASSUMED STRESS HYBRID FINITE ELEMENT FOR BUCKLING ANALYSIS

An assumed stress 8-node hybrid finite element is presented for buckling analysis. The formulation is based on Hellinger-Reissner variational principle. The element has 24 membrane degrees of freedom and 24 bending degrees of freedom. The stress stiffness matrix is derived by using the Green-Lagrange strains. The element performance is assessed by comparison with numerical examples.


INTRODUCTION
Buckling loads have detrimental influence on the performance, operation and maintenance of structures in field of engineering, i.e. in civil, mechanical, aerospace, automobile and marine structural engineering.Buckling analysis is carried out to calculate the loading value at which the structure buckles.Theoretical solutions can only be available on some simple structures under simple loading and with simple boundary conditions.Finite element method has been considered as the most effective approach in buckling analysis and used extensively.Fafard et.al. [1] presented a general two dimensional thin plate/shell theory for the study of elastic stability.M. Talbot and G. Dhatt [2] presented discrete Kirchhoff elements for shell analysis with large geometrical non-linearities and bifurcations.Rengarajan et.al. [3] developed a shell finite element for linear stress, buckling and free vibration analyses, Li et.al. [4] presented a shell element by a special transformation solid element for buckling analysis of thin-walled structures.Kim and Kim [5] developed practical design equations estimating the buckling strength of the cylindrical shell and tank subjected to axially compressive loads by using the ABAQUS, a commercial finite element analysis program.
The presented paper places emphasis on buckling analysis of structures by using finite element method.An assumed stress finite element based on Hellinger-Reissner functional is developed.The efficiency, accuracy and versatility of the proposed element are demonstrated by some numerical examples.

ELEMENT FORMULATION
Since Pian [6] first established the assumed stress finite element model and derived the corresponding element stiffness matrix, the hybrid stress model has been shown highly accurate, and easy to fulfill the compatibility condition of the finite element method.The element developed in this paper is based on Hellinger-Reissner variational principle.The Hellinger-Reissner functional can be written as where {} is the stress vector, [S] is the material stiffness matrix relating strains, {}, to stress ({}=[S]{}), [D] is the differential operator matrix corresponding to the linear strain-displacement relations ({}=[D]{u}) and V is the volume of structure.
The approximation for stress and displacements can now be incorporated in the functional.The stress field is described in the interior of the element as {}=[P]{} (2) and a compatible displacement field is described by {u}=[N]{q} (3) where [P] and [N] are matrices of stress and displacement interpolation functions, respectively, and {} and {q} are the unknown stress and nodal displacement parameters, respectively.Intra-element equilibrating stresses and compatible (boundary or intra-element) displacements are independently interpolated.Since stresses are independent from element to element, the stress parameters are eliminated at the element level and a conventional stiffness matrix results.This leaves only the nodal displacement parameters to be assembled into the global system of equations.Therefore, use of hybrid-stress versus assumed-displacement elements can be made transparent to general-purpose program users.Substituting the stress and displacement approximations Eq. ( 2), Eq. (3) in the functional Eq.( 1) where The form of Eq (5.a) and Eq (5.b) is directly amenable to numerical integration (i.e.Gauss quadrature).In principle, integration rules should be chosen which integrate all terms in [H] and [G] exactly.Arbitrary reduction of the integration order in [H] can lead to poor conditioning or singularity in [H] so that [H] -1 will not exist, and reduction in [G] can lead to spurious zero energy modes since [G] controls element stiffness rank.Now imposing stationary conditions on the functional with respect to the stress parameters {} gives Substitution of {} in Eq. ( 4), the functional reduces to where is recognized as a stiffness matrix.
The solution of the system yields the unknown nodal displacements {q}.After {q} is determined, element stresses or internal forces can be recovered by use of Eq. ( 6) and Eq. ( 2).Thus In modelling structures using displacement based or hybrid elements, body forces applied to the elements are replaced by equivalent nodal forces.With this replacement, the stiffness and stress matrices in element formulations need only to be considered for forces applied at the nodes.
The stress stiffness matrix is derived by using the Green-Lagrange strains which correspond to defining the strain of a line segment by the equation where ds and ds * are respectively the initial and final lengths of the line segment.
The nonlinear terms are retained only for the membrane strains.The membrane strains are written as where the linear part is given by and the nonlinear part is given by The bending strains are computed using the changes in curvature by The transverse shear strains are given by   The generalized Hellinger-Reissner functional including the nonlinear strains can be written as where   o  is the prescribed prebuckling stress state.Substituting the stress and displacement approximations Eq. (5.a), Eq. (5.b) in the functional and stress stiffness matrix is given by Here  N is obtained from shape functions   N by appropriate differentiation and ordering of terms.
The matrix    , which corresponds to the pre-buckling stress state, consists of membrane stress resultants which are evaluated from a linear static stress analysis.It is given by where

FINITE ELEMENT FORMULATION
The geometry and node numbering for the element is shown in Fig. 1.At each node of this proposed finite element the unknown displacements are translations u, v, w and rotations  x ,  y ,  z .
In practice, initial polynomials are usually assumed for the stresses after which the equilibrium equations are applied to these polynomials yielding relations between the 's and ultimately the final form of [P].The equilibrium equations which are applied to stress field polynomials yielding relations between the 's is given in Eq. (24).
The stress field should be selected in such a manner that no spurious zero-energy mode is produced.A spurious zero-energy mode is produced when the product of a selected stress term and the strains that are derived from the displacement functions produces zero strain energy under a particular deformational displacement field.The number of stress parameters, which is the number of columns in [P], must be at least equal to the number of degrees of freedom of the element less the number of degrees of freedom necessary to prevent rigid body motion Darılmaz and Kumbasar [7], Darılmaz [8].Spurious zero-energy modes generally occur for regular geometries such as rectangular planar elements and brick solid elements and disappear for irregular geometries, Pian and Chen [9].
The membrane part of the element has 23 degrees of freedom and the following equilibrating stress resultant field is considered for this part.
The bending part has also 27 degrees of freedom and the following equilibrating stress resultant field is considered for this part.
The Reissner-Mindlin plate theory is used for the bending part.The bending part of the element requires only C 0 continuity and takes into account the effects of transverse shear deformations by assuming constant transverse shear strains through the thickness of the plate.A static form factor k=5/6 is included in the transverse shear strain energy that accounts for a parabolic distribution of shear stress over a rectangular section.
In this study, since the stress field exactly satisfies the equilibrium equations given in Eq.( 24) the evaluation of [G] matrix is carried out along the boundary of elements.The components of stress resultant field along the boundary are given in Eq. ( 27) here [T  ] is the stress resultant transformation matrix.The displacement components along the boundary is also given in Eq. ( 29).
The translational and rotational components on the element boundary are defined with respect to nodal values as in Eq. (31).
are the shape functions in coordinates (,).It is worth mentioning here that the true nodal normal rotations are given by 1/2(v/x-u/y) evaluated at the nodes.Hence the terms  zi are not true nodal rotations.
Advantages of including rotational degrees of freedom in a membrane element are as follows.Major finite element programs allow six degrees of freedom per node.Accordingly, since rotational freedoms are already available in the program structure for frame elements.Including normal rotation as a nodal degree of freedom in flat shell elements is a neat way of avoiding the zero-stiffness condition that arises when all shell elements at a node happen to be coplanar.Finally, element performance may improve.
The differentiation of the work done by the distributed load with respect to the nodal displacements gives the load matrix The distributed load {F} can be expressed in terms of nodal load values {F} e as in Eq. ( 34) Substitution of Eq. (34) into Eq.( 33) , yields

Buckling load of a beam
The buckling load of the beam having a channel section is determined.The material properties are, elasticity modulus E=2x10 8 kN/m 2 , Poisson's ratio =0.3.Dimensional and geometrical properties of the beam are depicted in Fig. 3. Comparisons of the results of the ANSYS (shell93) and proposed element with the theoretical solution are shown in Table 1.From Table 1, the agreement is seen to be good.Buckling mode of the beam is given Fig. 4.  [11] 147.

Buckling of a square plate
The elastic stability of a square plate structure subjected to inplane edge load is analyzed.Material properties of the plate are taken as elasticity modulus E=3x10 7 kN/m 2 , Poisson's ratio =0.3.Its geometrical characteristics are shown in Fig.

Folded plate structure
The elastic stability of a folded plate structure subjected to distributed load is analyzed.Material properties of the plate are taken as elasticity modulus E=2x10 7 kN/m 2 , Poisson's ratio =0.2.Its geometrical and material characteristics are shown in Fig. 6.

Cylindrical Shell
The elastic stability of a cylindrical shell subjected to distributed external load is analyzed.Its geometric characteristics are shown in Fig. 8 with radius R= 1 m, height h=2 m, and thickness t = 5 mm.As a boundary condition bottom of cylinder is fixed and top of cylinder is free.It is made of steel with the elasticity modulus E =2.1x10 8 kN/m 2 and the Poisson's ratio =0.3.The computed elastic bifurcation buckling load obtained by using the presented element is compared with ANSYS solution and given in Table      [11] 8412.7 16X16 13 Like previous examples, for this example presented element solution is in a good agreement with reference solution.

CONCLUSIONS
The development of an assumed stress hybrid element is presented for buckling analysis.The element is developed through the combination of the membrane element with a drilling degree of freedom and the Reissner-Mindlin plate bending element.Thus, the combined element element possesses six d.o.f per node, which is very conveniently used in combination with other six d.o.f per node elements for giving more accurate solutions.On the basis of the representative numerical examples, the good accuracy of the proposed element for buckling analysis of structures has been demonstrated.From the numerical results of previous section, the accuracy of the element is satisfactory and can be implemented easily to analyze elastic stability problems.
any row of the matrix [G] are the forces at the nodes when the stress in the element is represented by one column of [P].The elements in any row of [G] represent a set of self-equilibrating forces.

Fig. 1
Fig.1 Geometry and node numberingThe stress field and the body forces have the following components

Fig. 2
Fig.2Components of stress resultants along the boundary of an element and defined as in Eq. (28) [T u ] is the displacement transformation matrix.Substitution of Eq.(28.a) and Eq.(29.a)into Eq.(5.b) , yields

Fig. 3
Fig.3 Dimensional and geometrical properties of the beam

Fig. 4
Fig.4 Buckling mode of the beam under axial load.

Fig. 5
Fig.5 Dimensional and geometrical properties of the plate .

Fig. 8
Fig.8 Cylindrical shell and its buckling mode under external load.During the analysis, 16 elements are used along the circumference and height of cylinder.

Table . 1
Elastic buckling load of the beam (kN)

Table . 2
Elastic buckling load of the square plate (kN) Solutions for first four buckling modes are compared with reference solutions and close agreements are obtained.Results are given in Table.2.

Table . 3
Fig.6 Folded plate structureSolution is compared with ANSYS solution and a close agreement is obtained.Results are given in Table.3.The first two buckling modes are given in Fig.7.Elastic buckling load factors of the folded plate structure Fig. 7 Buckling modes of folded plate structure

Table . 4
Elastic buckling load of a cylindrical shell ( kN)