The family of multilayered finite elements for the analysis of plates and shells of variable thickness

. Urban development requires careful attitude to environment on the one hand and protection of the population from the natural phenomena on the other. To solve these problems, various building structures are used, in which slabs and shells of variable thickness find the wide application. In this work, the family of multilayered finite elements for the analysis of plates and shells of variable thickness is described. The family is based on the simplest flat triangular element of the Kirchhoff type. The lateral displacements in this element are approximated by an incomplete cubic polynomial. Such an element is unsuitable for practical use, but on its basis, improved elements of triangular and quadrilateral shape are built. Particular attention is paid to taking into account the variability of the cross-section. The results of the developed elements testing are presented, and the advantages of their use in the practice of designing and calculating the structures are shown.


Introduction
Plates and shells of variable thickness are used in various fields of technology, both as the independent elements and as the part of combined systems.In aviation, for example, such structures include blades of propellers of airplanes having a complex spatial configuration, and some elements of the wing design [1][2][3].In the construction the slabs and shells of variable cross-section are used mainly in reinforced concrete structures.A classic example of a shell of variable thickness is, for example, a cooling tower shell [4].Variable crosssection panels are used in frameworks of bridges, retaining walls, overpasses, prefabricated large-size slabs for industrial and civil buildings of types "T", "2T" and other structures [5][6][7][8][9][10].In particular, when designing reinforced concrete floor slabs or roof slabs, a method of stress reduction is often used both in the marginal and in the main part of the slab due to its gradual thickening towards the contour line.In some cases, this technique significantly changes the working conditions of the plate, which begins to work like a depressed shell, as a result of its middle surface warping.When calculating such constructions using the finite element method, it is necessary to use shell finite elements of both constant and variable thickness.However, variable thickness elements are available not in all finite element programs.Therefore, calculators and designers using such programs are forced to simplify design schemes and present slabs of varying thickness as stepped slabs [11][12][13][14][15][16][17][18].This leads to a distortion of the stress-strain state of the structures being calculated.It should be noted that in the works devoted to the calculation of reinforced concrete slabs taking into account the non-linearity of deformation, plates of constant section are also mainly considered (see, for example, [19][20][21][22][23][24][25]).The lack of accessible and approved methods for calculating plates of variable thickness creates problems in the design and calculation of structures containing similar elements.It is possible that this is one of the reasons that, in practice, plates of constant thickness are used in cases when there is no need from the point of view of the logic of the structure work, for example, when arranging balconies or loggias, because it leads to waste of material.The above-mentioned problems are solved in the computer program PRINS, which includes finite elements of plates and shells of variable thickness of triangular and quadrilateral shape, intended for both linear and non-linear calculations.

Methods
The following describes the family of multilayer finite elements implemented in the computer program PRINS, a general description of which can be found in [22,23].This article discusses the features of the forming of stiffness matrices associated with the variability of the cross-section of elements.
The basis of the developed family is the simplest multilayered element of a triangular shape, shown in Fig. 1.The first degree polynomials for displacements in the plane and the incomplete cubic polynomial for lateral displacements are taken as approximating functions.It is assumed that the thickness of the element varies linearly.The division into layers, which was adopted during the formation of element characteristics, is shown in Fig. 1.The stiffness matrix of an element is calculated by the formula (see, for example, [16]); where Bis the geometric matrix, connecting the components of the deformations and displacements, andC is the physical matrix, connecting the stress components with the components of the strains.When using the Kirchhoff hypothesis and adopted approximating functions of displacements, the matrixB is reduced to the form [16]: where B P , A P −1 and A b −1 −some numerical matrices, B b − the matrix, whose elements depend on the x and y coordinates of the inner points of the finite element (methods for calculating of these matrices are well known; see, for example, [16]).Substituting (2) into (1), we obtain where We now take into account the variability of the thickness of the element.We consider the course of calculations on the example of the K b matrix.
Integration over the coordinate z in the second line of formula (4) will be performed in layers, assuming that the properties of the material do not change over the thickness of the layer.We get: (5) where  , and  , are the coordinates of the lower and upper boundary of the layer along the z axis.
The integration in the formula ( 5) is performed over the area F of the triangle.
Values  , and  , depend on the coordinates x and y.We assume that the thickness of the plate varies according to a linear law.Then for any layer of the slab, we can write: The coefficients  1 ÷  3 entering into formula (6) can be found from the condition that the z coordinates of the layer at the nodal points are known.We get: -at  = 0 and  = 0   =  , =  1 ; -at  =   and  = 0   =  , =  1 +  2   ; -at  =   and  =     =  , =  1 +  2   +  3   .
From the above conditions we find: With account of relations (7), formula (6) In accordance with (8), the z-coordinates of points of the upper and lower surfaces of the layer will be determined by the formulas: We write the formulas ( 9) and (10) in the form: The expression under the integral sign in the formula (5) represents the matrix of 9 × 9 order.We denote this matrix as follows: ) B b .
The expression for the matrix K b takes the following form: The calculation of the matrixK b in general form is difficult, so the integral ∫ R F dfis found by numerical methods.
The triangular element considered above has well-known drawbacks (see, for example, [16]).But on its basis improved elements of both triangular and quadrangular shape are obtained.
The improved triangular element is formed by averaging of the characteristics of the three sub-triangles described above.The local coordinate system is introduced for each of the subtriangles, as shown in Fig. 3.And initially all of the main characteristics of elements are calculated in this systems.The coordinate system of the first sub-triangle is taken as the local coordinate system of the composite triangle.Then all of the characteristics of sub-triangles are converted to the axes of the first one by standard transformations.For example, the stiffness matrix of a composite triangle is found by the formula where K 1 , K 2 and K 3 are the stiffness matrices of sub-triangles 1, 2 and 3 in their local axes; L 2 and L 3 -matrixes of the direction cosines of the axes x 2 y 2 and x 3 y 3 in the axesxy, respectively.
Similarly, other characteristics are found and averaged.
To obtain the stiffness matrix of a quadrilateral finite element, we use the technique of [24].We divide the quadrilateral (not necessarily rectangular) finite element into triangles, as shown in Fig. 5.The stiffness matrix of a quadrilateral element in its local axes is found by averaging of the stiffness matrices of the triangles reduced to these axes, according to the formula: The stiffness matrices of the triangles in their local axes are found by the formula (1), and the transition to the local axes of the quadrilateral is carried out by the standard transformation:

Results
To verify the accuracy of the developed family of elements, a number of test problems were solved.Two of them are listed below.Task 1. Cantilever plate of variable thickness, elongated in one direction, was calculated.The dimensions of the plate are shown in Fig. 6.Two variants of calculation schemes were used with a grid of nodes4 × 4 and4 × 12.
Figure 7 shows the stresses in the finite elements on the upper surface of the plate for the second variant.The calculation results are summarized in Table 1 (the stresses on the plate surface were averaged over a series of elements with the same "y" coordinate).The analytical solution is obtained by the strength of materials methods with a modulus of elasticity = 3,2 × 10 7 kPa.
The error in determining of the displacements with a grid of nodes4 × 4 was equal to 2,38%, and with the grid of nodes 4 × 12to 0.79%.The maximum error in the determination of stresses with a grid of nodes4 × 4 is equal to 9,8%, with a grid 4 × 12 -0.85%.The resulting error is quite acceptable for engineering calculations.Task 2. A plate of triangular shape, having dimensions and loaded, as shown in Fig. 8, was calculated.The material used was concrete with the following characteristics:  = 3,2 × 10 7 kPa; = 0,2.The plate was clamped on the left side.
The calculation scheme of the finite element method is shown in Fig. 9, and the stresses on the upper surface are shown in Fig. 10.Stresses were drawn at the center of gravity of each triangular element.To assess the accuracy of the developed element, an analytical calculation was made using the beam theory.The slab was considered as a beam of variable stiffness.Displacements were determined from the differential equation which was formulated in axes  1  1 (see fig. 8).The moment of inertia of the cross-section was determined by the formula ( 1 ) 1 = ( 1 )ℎ 3 ( 1 ) 12 .
At the given dimensions of the slab, the maximum displacement found from equation (a) was equal to 0.2637 m.The maximum displacement in the finite element analysis of plate is 0.267 m.The difference in the maximum displacement was 1.25%.
The analytical values of the stresses were assumed to be uniformly distributed over the width of the section and were calculated by the formula The average stress value  in the centers of gravity of a number of elements adjacent to the embedding, obtained by finite element method, was 14766 kPa.The value of 14461 kPa was found by the formula (b).The discrepancy was 2.1%.The stress at the center of gravity of the element adjacent to the wedge apex, calculated by finite element method, is 23200 kPa.The value  , kPa was found in this section by the formula (b).The difference in stresses in this section is 0.93%.
Thus, the developed element provides acceptable accuracy of the results both in displacements and in stresses.As an illustration of the capabilities of the program PRINS in the design of structures, we consider the calculation of a cantilever plate for a uniformly distributed load and its own weight.
The plate was calculated in two versions.In the first variant, the plate thickness was taken constant, in the second -variable.The dimensions of the plates are shown in Fig. 11 and 12.
The calculation was carried out in a linear formulation with the following input data: heavy concrete of class B40 with the initial module = 33200 MPa, Poisson's ratio = 0.2, weight density  = 25 kN/ 3 and tensile ultimate strength   = 1,4 MPa.The plates were reinforced in the stretched zone by A400 class reinforcement with a cross-sectional area equal to 62,8 sm 2 /.The intensity of the load was taken equal to 1 кPa.
For both options, the same calculation scheme was used (see Fig. 13).The results of the calculations are presented in Fig. 14 -15.

Discussion
As can be seen from Fig. 14 and Fig. 15, both calculated plates give almost the same stress values on the stretched surface in the rigid support:  kPa х, in the case of a plate of constant cross-section and  kPa х, in the case of a plate of variable cross-section.Therefore, the compared plates provide the same strength.However, the volume of a slab of variable cross-section is much smaller.With the accepted dimensions, the volume of the slab of constant cross-section is   = 1,2 m 3 , and the volume of the slab of variable crosssection is  m 3  .Thus, the replacement of a slab of constant cross-section with a slab of variable cross-section provides in this case more than twice the savings of concrete.
As for deformations of the plates, the maximum deflection in the case of a plate of constant cross-section was equal to  ,  = 0,511 mm, and in the case of a plate of variable cross-section  mm , , i.e. a bit more.However, in the latter case, the maximum deflection is 1/2381 of cantilever length, equal to 3 m, which is much less than the norms adopted in construction.
On the basis of the calculations performed, it can be concluded that the program PRINS may be useful in the design and calculation of structures containing plates of variable thickness.

2 , 3 , 4 where
K m xy −the stiffness matrix of the m-th triangle in the local axes of the quadrilateral, K m − the stiffness matrix of the m-th triangle in its local axes, L m − the matrix of the direction cosines of the local axes of the m-th triangle in the local axes of the quadrilateral.

Fig. 7 .
Fig. 7. Stresses on the top surface with a grid.

Fig. 11 . 12 .
Fig. 11.The plate of the constant Fig. 12.The plate of the variable thickness.

Fig. 14 .
Fig. 14.Stresses on the upper surface of the plate of constant cross-section (kPa).

Fig. 15 .
Fig. 15.Stresses on the upper surface of a plate of variable cross-section (kPa).

Table 1 .
Comparative analysis of the results.