Base functions effect in stability finite element analysis of compressed plate

Geometric nonlinear solution of a compressed plate is presented in this paper. Basic assumptions are specified and incremental conditional equations are derived from the variational principle of minimum of total potential energy. Full Newton-Raphson procedure, in which the stiffness matrix is updated at every equilibrium iteration, has been applied for solution. The importance of modifications of base functions for solving geometric nonlinear problems is analysed. The solved example is presented, the differences are compared and explained.


Introduction
Thin rectangular plates are the simplest shapes, easy to made therefore they are widely used in many branches of engineering. They occur as parts of aircraft, marine and vehicle structures in mechanical engineering, form bearing parts as well as reinforcements in civil engineering.
Thin plates subjected to the external pressure are liable to the buckling due to dominant compression membrane forces within the plate. Stored strain energy can be suddenly released and converted into kinetic energy with subsequent collapse of the structure. It is the reason why the stability problem has been analyzed since the beginning of the twenty century. Solving postbuckling behaviour of the compressed plate, nonlinear analysis is essential. It is also necessary to include initial imperfections of a real plate into the solution. Experience shows that small deviations from the prescribed geometry of a perfect plate may cause a strong reduction of the load-carrying capacity.
In this paper the effect of choice of base functions upon the mesh density of thin plate structures as well as upon the convergence of FEM solution of stability problems, have been studied. Numerical analysis of stability problems has been converging from above, therefore the mesh density should be identified at which the obtained results are applicable. Nonlinear solutions, to the contrary (incremental, refined by iterations) are time-consuming and the dense meshing is frequently inapplicable in practice. The paper is focused on a search for a proper type of a finite element. By sequential adding the degrees of freedom one can achieve higher accuracy using smaller number of iterations. Our aim is to find an optimal ratio between the speed of solution convergence and the degree of difficulty of the calculation, to find the proper element for stability problems of thin plates, applicable also at small number of finite elements.

Theoretical assumptions
Geometrically nonlinear theory deals with more complicated strain -displacement equations. Let us assume a rectangular thin plate with the thickness t simply supported along the edges (Fig. 1). The displacements of the point of the neutral surface are denoted q = {u, v, w} T and the related load vector is denoted p = {px, 0, 0} T . Assuming the so called von Kármán theory, the normal displacements (w) are much bigger than tangential displacements (u, v). By formulation of the strains, it has to be taken into account nonlinear terms.
, the indexes denote the partial derivatives.

Configuration of the System of Conditional Equations
The increments of displacements of the point of neutral surface are denoted and the increments of the strains can be written as: Restricting to the isotropic elastic material and to the constant distribution of the residual stresses over the thickness, the increments of the stresses can be expressed as The initial displacements will be assumed as the out of plane displacements only. For such defined task, one can express the increment of the total potential energy as: The system of conditional equations one can get from the condition of the minimum of the increment of the total potential energy 0 U =   . After modifications this system can be written as: where is the incremental stiffness matrix of the plate,

Solution of the system of conditional equations
If we used a rectangular element obtained by composition of eight in-plane parameters and twelve plate parameters, we would have to choose high density of meshing to finite elements, in order to receive sufficiently suitable (accurate) solution. In other words, the solution would converge very slowly. In many papers [e.g. 2,3] this phenomenon was demonstrated and illustrative examples are presented. Hence, for stability problems of thin plates was used a modified rectangular four-node finite element with 48 degrees of freedom (12 parameters in each node) with following choice of base functions [4]. The solution of the system (5) cannot be achieved directly, and is based on a step by step incremental process that causes a deviation from the equilibrium nonlinear path. To correct this deviation, an iterative technique, Newton-Raphson method is used [5]. In the case of a more complicated loading path, it is possible to utilize a calculation of controlled displacement increment [6]. In the case of snap-back, this option fails, but there are still some perspectives in the possibility to change the pivot parameter during the calculation. In this case the interactive approach and calculation check are requisite.

Design of computational model and numerical simulation
Numerical analysis of steel (E = 210 GPa, ν = 0.3) rectangular plate a = b = 500 mm, performed by classic NRM controlled by load increment. Constant value Δp = 5 kN/m was considered for all types of finite elements and for each mesh density. The Ansys system [7] was used for analysis, these elements were stepwise used to create the model: A) Shell 181 -classic bilinear element, B) Shell 281 -8 nodes biquadratic element with midside nodes, C) User 33 -element with 48 D.O.F. according to (5) -aproximation with bicubic polynomial. In all cases the plate was modeled with 2x2, 4x4, 8x8 and 16x16 elements respectively (see Fig. 2).   If we compare a number of equations that must be created and solved in each iteration (which also represents a dimension of matrix Kinc, which has to be repeatedly factorized) we realize that: while in cases A and B the numbers are comparable, dimension C is smaller and the solution is faster. Utilizing of this element in the solution therefore appears to be more favorable. To formulate a more general conclusion, one can say that making the solution more accurate by increasing the number of nodes in element is not a very appropriate option. Increasing of the number of nodal parameters seems to be preferable for this purpose.
Presented results have been arranged due to the research supported by the Slovak Scientific Grant Agency, project No. 1/0265/16.