DEVELOPMENT OF A GENERAL ALGORITHM FOR SOLVING THE STABILITY PROBLEM OF ANISOTROPIC PLATES

, to solve stability problems of anisotropic plates, different representations of the bending deflection function in different rows are used. But the use of such representations is justified only under certain boundary conditions and under the condition of uniformly distributed load. The study described in this paper offers a way to overcome these difficulties, allowing the numerical values of critical forces to be determined without much difficulty. With increasing grid density, the accuracy of the critical load value increases rapidly and with an 8 × 8 grid, the deviation from the exact solution equal to is 1 %. From a practical point of view, the discovered mechanism of numerical realization of this problem allows to improve engineering design calculations of stability of anisotropic plates with different conditions on supports and with different loading


Introduction
The improvement of engineering structures, the complication of their operating conditions, as well as the introduction of new materials lead to the need for further development of methods for calculating structures.With the use of new structures and materials, material consumption is reduced, the scope of application of thin-walled systems with low rigidity, for which the danger of elastic loss of stability increases, and, consequently, the importance and relevance of the theory and methods for practical solution of problems of elastic stability of engineering structures increases.The problem of stability is relevant to all types of structural elements and machine parts, and its importance especially increases in relation to anisotropic thin plates.The practical solution of this problem currently encounters difficulty mainly not in drawing up equations, but in obtaining the numerical value of the critical parameters of the external influence.
The approximate method in general and, in particular, the discrete method for solving the problem of stability of elastic systems has a certain area of application, the main indicators of which follow from the fundamental characteristics of the method.First of all, these include the necessary initial information.In this method it is an additional load.To obtain it, one must have differential equations that describe the deformation of the system in a given and adjacent states.Since such a condition meets almost no restrictions when setting stability in the "small", the scope of application of the additional load method is very extensive.
Therefore, research on the development solution the problem of stability and the discrete method for solving is relevant for all types of structural elements and machine parts.

Literature review and problem statement
The works [1,2] of many authors, classics of the mechanics of deformable solids, are devoted to the stability of elastic systems, which analytically consider stability examples following which it is difficult to determine numerical values

This paper is devoted to the development of general algorithm for solving to the stability problem of anisotropic plates using the additional load discretization method. The study of the stability problem is relevant for all types of structural elements and machine parts, and its impor-
tance is especially increasing with respect to anisotropic thin plates.This is due to the fact that with the use of new structures and materials, the material intensity is reduced, the area of application of thin-walled systems with low stiffness, for which the danger of elastic loss of stability increases, and, therefore, the importance and relevance of the theory and methods of practical solution of problems of elastic stability of such structures increases.

In many works, analytical expressions for determination of critical load are given. At present, the determination of critical loads causes great difficulties in their numerical determination. Therefore, the article presents the most effective numerical and analytical solution of this problem. As a rule, to solve stability problems of anisotropic plates, different representations of the bending deflection function in different rows are used. But the use of such representations is justified only under certain boundary conditions and under the condition of uniformly distributed load. The study described in this paper offers a way to overcome these difficulties, allowing the numerical values of critical forces to be determined without much difficulty. With increasing grid density, the accuracy of the critical load value increases rapidly and with an 8×8 grid, the deviation from the exact solution equal to is 1 %. From a practical point of view, the discovered mechanism of numerical realization of this problem allows to improve engineering design calculations of stability of anisotropic plates with different conditions on supports and with different loading Keywords: structural stability, plate stability, plate strength, anisotropic plates, orthotropic plates, discretization method
of the critical load.The paper [3] considered the solution of the stability problem for shallow shells by the method of discretization the additional load and the basis of the method of discretization the additional load in relation to the problem of stability of curved rods of arbitrary shape.A numerical solution of the stability problem is shown, but only problems related to curved rods are solved here.Some issues directly related to the method of discretization the additional load are considered by the author [4,5] in relation to the current state of the question of the stability of elastic systems.The numerical solution is also shown here, but the anisotropic material characteristics of the structure are not considered.The paper [6], have considered the problems of calculating the stability of rods and plates depending on their geometric parameters and analyzed the obtained calculation results, determined the factors influencing the stability of rods and plates and the limits of applicability of formulas for calculating the critical buckling force of rods and plates.However, the paper does not show the complex stability problems of plates with different boundary conditions.
The paper [7] shows a method for calculating the stability of plates and rods, which allows for various calculations of the stress-strain state depending on the boundary conditions -applied loads, moments and degrees of freedom.The work is interesting from an engineering point of view, but the non-uniformity of the distributed compressive load is not taken into account.
In article [8], the calculation of the stability of an isotropic plate by the method of separation of variables was considered, the plate was calculated for stability under compression in the longitudinal direction, and the results were compared with existing solutions.The variable separation method was used, however, the paper does not consider the anisotropy of plates often encountered in engineering structures.
A method for calculating the stability of plates under shear loads is presented [9], but the paper does not give numerical solutions of these equations for non-uniform loading.Based on Vlasov's variational method, a system of differential equations was obtained to study the stability of plates.In [10], the problem of stability of a compressed plate is solved and the theoretical results are compared with data from a full-scale experiment, but the paper does not take into account the anisotropy of the plates.A method has been proposed [11] for determining the critical load leading to loss of stability of anisotropic rectangular simply supported plates, complex plate supports are also not considered here.
In the reviewed papers [1][2][3][4][5][6][7][8][9][10][11], many stability problems of structures have been solved, but the issues related to the numerical implementation of the problem of determining the critical force have not yet been solved.The reasons for this may be objective difficulties related to the complexity of solving stability problems.
It is clear that a general method for numerical solution of the stability problem of anisotropic plates with a variety of boundary conditions and contours under different loads is needed.

The aim and objectives of the study
The aim of the study is to develop a general algorithm for solving the stability problem of anisotropic plates with different boundary conditions and different external load.This will allow to solve in practice the stability problems of anisotropic plates as a part of various structures.
To achieve this aim, the following objectives are accomplished: -to develop a solution to the problem of stability of anisotropic plates using the method of discretisation of additional load, -to determine the Green's function by the variational-difference method as an invariant part of the problem of determining the critical load; -to show the convergence of the numerical implementation of this problem using the proposed method of additional load discretisation.

Materials and methods
The study deals with the development of a solution to the problem of stability of anisotropic plates.The elastic stability of plates is investigated.The material of the plate is assumed to be elastic.The hypothesis of straight non-deformable normals is also used.Small summands of the second order of insignificance are neglected in the differential dependences of deformations and relative deformations.
The stability of a deformable system is considered as a special type of its ultimate resistance, which can occur before the system exhausts mechanical resistance according to the strength conditions of the material.The design of structures for stability consists of determining the critical load corresponding to this ultimate resistance.It determines the load-bearing capacity of a structure based on stability conditions.During the operational stage, the load value must be below the critical value in compliance with the normalized stability margin.The task of the theory of stability of elastic systems is to formulate mathematically the stability condition, create and solve the resulting equation to determine the value of the critical load.
Determining the ultimate resistance as the final stage of the deformation process requires its broadest description.It can be obtained by nonlinear equations that express the deformation process more completely than a linear one.
In a stability problem, of the two nonlinearities -physical and geometric -the second one is mandatory, since the resistance limit under the stability condition precedes the resistance limit under the material strength condition.The basis of geometric nonlinearity is the theory of nonlinear deformations.
Let's consider the qualitative features of the linear deformation components, which are of significant importance for the problem of stability of deformable systems.In the nonlinear, as in the linear theory of deformations, two types of deformations are distinguished.The first one describes linear deformation: the bringing together, or removal, of the ends of an infinitesimal segment ds emanating from the point M(x,y,z) by a unit of its original length.It is determined in the same way as in linear theory by the formula: or by the formula: where ds 1 is the distance between the end points of the segment after deformation.
The second describes the shear deformation: the deformation of the angle between two infinitesimal segments intersecting at point M. It is also defined similarly to the shear strain in the linear theory of deformations.The peculiarity of nonlinear deformations, which distinguishes them from linear ones, stems from the fact that when expressing the deformation length of a segment through its corresponding length before deformation, rotations of the segment are taken into account, i.e. the change in the direction cosines of the segment ds is taken into account with its transition to ds 1 , which in the x, y, z axes are equal for the segment ds: l=dx/ds; m=dy/ds; n=dz/ds.For the segment ds 1 :l 1 =dx 1 /ds 1 ; m 1 =dy 1 /ds 1 ; n 1 =dz 1 /ds 1 Where dx, dy, dz and dx 1 , dy 1 , dz 1 are projections on the x, y, z axes of the segments ds and ds 1 , respectively.
Linear deformation e is a function of f(l, m, n) and is represented as: ( ) In particular cases, for segments dx, dy, dz parallel to the axes x(l=1,m=0,n=0), y(m=1,l=0,n=0), z(n=1,l= 0,m=0), from (4) follows the value of the coordinate components e, respectively: z zz Where approximately (using the Taylor series): After expanding the radicals ( 4)-( 6) into the Maclaurin series in ε ii , preserving only the first degree, let's obtain simplified approximate formulas, which are usually used in practical calculations: e x =ε xx , e y =ε yy , e z =ε zz .With the same approximation, the coordinate components of the shear deformation in the xy, xz, yz planes are determined by the formulas, respectively: , .
In ( 7)-( 12), the terms outside the brackets determine the deformations in the linear formulation; let's call them basic, and the terms enclosed in curly brackets define nonlinear additions to them; let's call them "additionals".Additionals necessarily accompany the main deformations, since through them the rotations of elementary segments emanating from the point M that naturally occur in a continuous medium are reflected due to any of the components of its movement u, v, w.
In this regard, "additionals" bring not only quantitative, but also (which is no less important) qualitative changes in that they: a) supplement the main deformation with new types of deformations; b) oblige (through rotations of segments) to take into account the deformed scheme of the environment, which is not taken into account in a geometrically linear formulation.
Such an account radically changes many ideas established by logic, based on the linear theory of deformation.For example, a rod centrally compressed along a straight axis, from the position of linear theory, experiences only uniform compression, and from a nonlinear position, compression is accompanied by bending, which is characterized by rotations of the axis of the rod: .
This kind of rotation significantly distinguishes the continuous medium, which is used in the mechanics of a deformable solid as a model of a structural material, from the regions of figures representing it in the understanding of Euclidean geometry.In the first, unlike the second, the deformation of any infinitesimal segment occurs not as an isolated, separate one, but, on the contrary, together with all the segments surrounding it.The deformed state that arises at a point is mathematically expressed by the strain tensor (which has complete analogy with the stress tensor).Therefore, any deformation of a segment must be considered in connection with the general deformation of the medium.Then it becomes natural to expect a change in the distance between the end points of a segment along the line of its original position in the general case not only due to the distance (as a "Euclidean" segment), but also due to rotation.
The equations of deformation of the medium in a geometrically nonlinear formulation are compiled in the same way as in a geometrically linear one -with the involvement of either static-kinematic conditions (namely, the conditions of equilibrium, continuity and physical law), or energy conditions, which are reflected in the variational principles (these include the principles possible movements and possible stresses).
In order to simplify the equations, an approximate account of the "additional" is allowed.It usually consists in the fact that some of the terms in the "additional" are omitted and its influence is not taken into account in all equations, and especially in the equations of the physical law.

Results solution to the problem of stability of anisotropic
plates using the additional load discretization method

1. Development a solving the stability problem of anisotropic plates by the method of discretization of additional load
To obtain the additional load, one must have differential equations that describe the deformation of the system in a given and adjacent states.Since such a condition encounters almost no restrictions when setting stability in the "small", the scope of application of the additional load method is very extensive.It includes not only simple elastic systems of individual elements, but also complex ones formed by the articulation of elements.It is important to note that only differential equations are needed, not their solutions.
The former, unlike the latter, are preserved under any continuous and discontinuous boundary conditions, including for holes.In this regard, in relation to multi-element systems, it is enough to have differential equations for each element separately even before the conditions for the junction between them are met.
To determine the Green's function in complex cases, approximate methods are used: variational-difference method, finite element method and other methods.They solve a wide class of plane and spatial problems.Moreover, the Green's function, by its nature, as relating to a unit load, is an invariant part of the calculation; it does not depend on the load of the system and is preserved for any loads specified on it.By varying the grid of points when determining the Green's function, it is possible to identify a pattern, the use of which greatly simplifies computational operations.
It follows from the fact that the boundary conditions under the action of a single concentrated force have a local influence.Therefore, the matrix of equations compiled using the variational-difference method for determining the Green's function is invariant with respect to the zone of internal points remote from the boundary.This pattern allows to compile standard matrices of coefficients of equations for the Green's function in relation to each given system.
In this paper one considers the stability of an orthotropic (as one type of anisotropy) plate under the conditions of a plane problem, which is compressed by forces N x , N y in the direction of the x, y axes and sheared by forces N xy in its plane (Fig. 1).
The tension is determined for an orthotropic plate by the formulas: Poisson's ratios μ x and μ y are related by the relation: .
Bending stiffnesses D x and D y and torsional stiffness D k are calculated using the formulas: Bending and torque moments and the transverse forces arising in the cross section are determined through the curvatures k x , k y and k xy : ) , In formulas (25), (26) D pr refers to the reduced cylindrical stiffness: ( ) For an orthotropic plate, the equilibrium equation taking into account the deformed state has the following form: The right side of equation ( 28) represents an additional load q*(x,y), additionally obtained by a plate compressed in two directions by distributed forces (N x , N y ) and shifted in its plane by a force (N xy ) in an additional deflected state during loss of stability: ( ) Following the additional load sampling method, a grid is applied to the plate, for example, 5×5, 10×10 or another, formed by lines parallel to the x, y axes.Let's replace the additional load within the grid cells with dimensions Δx, Δy with averaged concentrated forces.Let's write down expressions for deflections at the nodal points of the cell.To do this, let's replace integration in (30) by summation and, using unit displacements δ ij , obtain the following equations for deflections w i : , The derivatives in the load expression q*(x i , y i ) are written in finite differences, the order of which is not higher than the second.Having performed the operations belonging to (31), let's obtain a system of linear homogeneous equations for the deflections w i .Let's equate the determinant of this system to zero: where a ii are known numbers; λ is parameter to be determined.
Having expanded the determinant (32), it is possible to arrive at an algebraic equation for λ under simple loading, when N x , N y , N xy depend on one quantity.The smallest value of which gives the numerical value of the critical force.

2. Determination of the Green's function by the variational-difference method
The Green's function ( ) , , , , k k w x y x y with the help of which it is possible to determine the elements of the compliance matrix δ ij in (31), is found using the variational-difference method.According to Lagrange's variational principle, the variation of the total potential energy is zero: Variation of the work of an external unit force 1 P = applied at point k: Variation of the work of internal forces throughout the entire volume of the plate: .
Substituting ( 35) and ( 34) into (33): General formula for strain energy for plates: ( ) The deformation energy for anisotropic plates is obtained from the general formula (37): here E x , E y are the elastic modulus of the plate material under tension and compression along the x, y axes; G is a shear modulus of elasticity; µ x , µ y is the Poisson's ratios.
To obtain the formula for potential energy in deflections, deformation in formula (38) is expressed in terms of deflections w, using the dependencies arising from the hypothesis of direct normals of the middle plane: In this case: Let's express derivatives through finite differences and replace integration by summing the values Π i ΔxΔy, where the discrete values of Π i at points i.Then, according to (42), it is possible to obtain Π as a function of deflections w i .Let's substitute Π in (36) and group the terms with the displacement variations δw 1 , δw 2 ,…, δw n taken out of brackets.Let's assume that a unit force when determining the Green's function is applied at point k=1, then according to (36): Each expression in parentheses must be equal to zero, since equation (43) must be satisfied for any value.Thus, equation (43) turns into a system of linear algebraic equations, the number of which is equal to the number of required displacements w i : Having solved the system of equations ( 44), let's obtain a diagram of deflections from a unit force applied at point 1.Due to the symmetry of the Green's function, it is an influence surface that gives us the values δ 1j ( j=1,2,…,n).Next, let's apply a unit force to the plate at point 2. In this case, system (44) will not change, only the unit on the right side will appear in the second row.Having solved this system, let's obtain the surface of influence of deflections at point 2 and, therefore, δ 2j .In the same way, let's determine the surface of influence of deflections for all other points of the plate.

3. Numerical solution of the plate stability problem using the additional load discretization method
It is required to determine the critical value of the force for a longitudinally compressed uniformly distributed force N x , simply supported square plate (Fig. 2) at, N y =0, N xy =0, N x =const, μ=0.3.
Let's first determine the Green's function using the variational-difference method.Let's apply a grid with 4×4 cells on the plate (Fig. 2).It is possible to determine the specific energy of transverse bending of the plate Π k (k=1, 2, ..., 9) at each nodal point of the grid: and so on.Total energy Π: Here, points 1-9 are points lying inside the plate at the nodes of the division grid; points: 10, 14, 18 and 22 are points located at the corners of the plate; the remaining 12 points: 11, 12, ..., 25 are intermediate points, three on each of the four sides of the plate.The coefficients at Π express the value of the part of the grid area on the plate from which the strain energy is collected in determining the total energy.
Following the Lagrange principle, let's vary Π with respect to w k and draw up equations to determine the deflections w k when applying a unit force sequentially at points k=1, 2, ..., 9.When applying a unit force at the point k=1, let's obtain a system of linear equations (49), where the unknowns are the deflections at points k=1, 2, ..., 9: Having solved this system, let's obtain discrete values of the Green's function as a function of influence for displacements 1 w at point 1.If on the right side the free term of the system of equations 2 16a D enters the second line, then it is possible to obtain discrete values of the Green's function as a function of influence for displacements 2 w at point 2 and etc.The matrix of coefficients in equations (49) remains unchanged for all positions of the unit force.The displacements k w obtained as a result of solving equations similar to (49), when applying a unit force sequentially at points k=1, 2, ..., 9, are equal to the elements of the compliance matrix of the system of equations (31) in lines1, 2, ..., 9, respectively.
Having , k w and following (31), it is possible to compose a system of homogeneous linear equations for w k .From the condition that the determinant of this system is equal to zero, an algebraic equation for N x,kr follows, solving which let's obtain: The considered problem was solved with different division grids.The results shown in Table 1.
Considering the numerical values of the critical load for different division grids (Table 1), let's draw the following conclusions: 1.With increasing grid density, the accuracy of the critical load value increases rapidly and with an 8×8 grid, the deviation from the exact solution equal to 2. The mesh density for the plate when discretization the Green's function and discretization the additional load can be different.This makes it possible to reduce the order of the matrices without compromising the accuracy of the solution for the critical load; the Green's function is invariant, i.e.Once defined, it is used unchanged under various system loads.

Discussion of results development solution the problem of stability of anisotropic plates using the additional load discretization method
The advantages of this research in comparison with similar known ones are as follows.Usually, to solve stability problems of anisotropic plates, different representations of the deflection function at loss of stability in different rows are used.But the application of such representations is justified only under certain boundary conditions and under the condition of uniformly distributed load.When attempting to go beyond these limitations, the numerical realization of the problem of determining the critical load becomes much more complicated.Objective difficulties arise due to the uncertainty of the mechanisms for representing the deflection function, taking into account the boundary conditions and the non-uniformity of the distributed load.The research described in this paper proposes a way to overcome these difficulties.It is based on the fact that the procedure of deflections function representation and boundary conditions consideration is carried out in discrete form, Green's functions (influence functions) are used to determine the additional displacements arising at the loss of stability.And Green's functions, due to their invariance, do not change when the external load changes.This method allows to obtain stable results for anisotropic plates with very different boundary conditions and outlines under different loading.
In contrast to the studies given in [4,5,10], the stability problem of an orthotropic plate, as one of the types of anisotropy, is solved on the basis of the same data as the stability problem of isotropic plates.
The solution of the stability problem by the method of discretization of additional load has shown that the proposed method, in contrast to the study [11], allows to determine the critical load by a general algorithm for arbitrary anisotropic plates, without complicating the solution of the problem for the most complicated cases of loading and boundary conditions.
The solved problems show the convergence of the numerical implementation of this problem using the proposed method of discretization of the additional load.The obtained value of the critical force is compared with the solutions of other authors.The comparison gives satisfactory results (Table 1).
In solving the stability problem of anisotropic plates by the proposed numerical and analytical method of discretization of the additional load, the determination of the Green's function (influence function) is of great importance.The Green's function, inherently related to the unit load, is an invariant part of the calculation; it does not depend on the load on the system and is preserved at any load given on it (42), (44).By varying the grid of points in determining the Green's function, a regularity has been identified, the use of which greatly facilitates computational operations.It follows from the fact that the boundary conditions under the action of a single concentrated force have a local influence.Therefore, the matrix of equations (49) compiled by the variational-difference method for determining the Green's function is invariant with respect to the zone of internal points distant from the boundary.This regularity allows to compile standard matrices of equation coefficients for the Green's function with respect to each specific system.
The peculiarity of the proposed method and the results obtained in comparison with existing methods [6,7,9] is that this method can be applied without difficulty to anisotropic plates with a large variety of boundary conditions and contours under different loading conditions.
The limitations inherent in this study are the limits of applicability of the additional load discretization method, such as solving the dynamic stability problem [7].
The disadvantage of this research is that increasing the accuracy of the solution leads to computational difficulties.
The development of this study can be in the field of stability of plates with different holes, stiffening ribs, in contact with an elastic base, etc.In these cases, the difficulty will be to correctly account for these variations.

Conclusion
1. Solving the stability problem by the method of discretization of the additional load has shown that this will allow to determine the critical load for arbitrary anisotropic plates without complicating the solution of the problem for the most complicated cases of loading and boundary conditions.
2. The Green's function determined in the study by the variational-difference method is an invariant part of the problem of determining the critical load; it is independent of the load on the anisotropic plate and is preserved under any loads given on it, which is one of the advantages of the proposed additional load discretization method.
where N x , N y , N xy are the forces in the plane of the plate in a given state.The deflection w(x k , y k ) at point k from the additional load is found using the Green's function(  )    , , , k k w x y x y : Fig. 1.Plate loaded in its plane

Table 1
Critical load values for different division grids