The New Approach to Analysis of Thin Isotropic Symmetrical Plates

: A new approach to solve plate constructions using combined analytical and numerical methods has been developed in this paper. It is based on an exact solution of an equilibrium equation. The proposed mathematical model is implemented as a computer program in which known analytical formulae are rewritten as wrapper functions of two arguments. Partial derivatives are calculated using automatic differentiation. A solution of a system of linear equations is substituted to these functions and evaluated using the Einstein summation convention. The calculated results are presented and compared to other analytical and numerical ones. The boundary conditions are satisﬁed with high accuracy. The effectiveness of the present method is illustrated by examples of rectangular plates. The model can be extended with the ability to solve plates of any shape.


Introduction
Thin plates are widely used as elements of building and engineering structures. The Kirchhoff plate theory is the basis of engineering design and calculation of such elements. It is an approximate zeroth-order shear deformation plate theory. Discrepancy between order of basic differential equation and number of boundary conditions is the main disadvantage of the Kirchhoff theory. There are considered only three from six physical equations. Shearing and normal transverse stresses are not taken into account. They are determined from equation of equilibrium. The higher-order shear deformation plate theories are free from this discrepancy.
According to the Kirchhoff plate theory solution of thin isotropic plates is equivalent to solution of fourth-order partial differential equation [1,2] for given boundary conditions. In the above equation w is the deflection of the plate, and D is its flexural rigidity. Function q(x 1 , x 2 ) describes distribution of load on upper surface of the plate. Uniform equation corresponds to (1) has the form In many cases, Kirchhoff theory gives satisfactory results and sufficient accuracy for practical purpose and correctly describes behavior of the structure. It is based on the linear distribution of the displacements through the plate thickness. As a result, parabolic distributions of the shearing stresses are obtained. For comparison Mindlin plate theory [3] also takes the same distribution of displacements but gives incorrect uniform distribution of the shearing stresses through the plate thickness. Higher-order theories like Christensen, Mindlin, Timoshenko can represent the kinematics better. However, they involve higher-order stress resultants that are difficult to interpret physically. Therefore, such theories should be used only when necessary [4]. In addition to its inherent simplicity and low computational cost, the first-order plate theory often provides a sufficiently accurate description of the global response for thin to moderately thick plates. Higher-order theories provide a small increase in accuracy relative to the first-order plate theory solution, at the expense of a significant increase in computational effort [5]. Therefore, the choice of the Kirchhoff theory to analysis of thin-walled engineering structures is still justified.
Let us note that exact solution of boundary problems of theory of the plates is possible only for simplest cases. At present, there are many various methods to solve thin plates, which can be divided into two basic groups: analytical [2,[6][7][8] and numerical ones [6,[9][10][11].
Analytical methods are more exact than numerical methods but they allow the solving of boundary problems of plates with a small number of unknown parameters and constrained by canonical contour as rectangle, square, triangle, circle or ellipse. Numerical methods are less exact but they allow the solving of arbitrary configuration plates [12] at the expense of a significant increasing of unknown parameters. It increases the time for calculations and leads to the accumulation of computational errors.
Numerical methods are currently reaching their limits. Among numerical methods, the Finite Element Method (FEM) [13,14] and Finite Boundary Method (FBM) [15][16][17][18] are the most common. The finite element method allows us to approximately solve plate of arbitrary configuration, where the analytical solution is impossible. Research on this method can be conventionally divided into: theoretical (related the construction of finite element model) and applied (related to the use of the method to solve complicated technical problems).
In recent years, many methods which are free from disadvantages of standard numerical ones have been suggested. They create separate group called analytical-numerical methods. According to these methods part of equations are performed in analytical manner and others with help of numerical procedures. The method suggested in present paper belongs to this group. Differences compared to numerical methods:

•
Simplicity of structure modeling. Equations of equilibrium are solved once independently of the boundary conditions and solution remains unchanged throughout the calculation process. This allows the reduction of structure modeling to set kinematic and static boundary conditions in the separate points at the plate contour. • Possibility to provide better accuracy in comparison with finite element method. It is achieved only by increasing of number of the approximations (parameter K) or replace nodes at the contour. This procedure is performed automatically with author's program to solve of the plate structure. • Possibility to define the load as a function of two variables. Because not all FEA programs provide such functionality, ABAQUS software have been chosen to compare the results.

•
Possibility of implementing of suggested method in higher-order theory, say Mindlin-Reissner theory.
The mathematical model proposed in the article has been implemented as a computer program which is characterized by a readable and concise source code (Appendix B). This method requires much less computation operations to solve the structure in comparison to FEM. There is no need to: • Discretize of the structure into set of finite elements, • Aggregate of this set into discrete structure, • Introduce local and global coordinate system, • Numerate of nodes. Solution of equilibrium equation can be obtained with given accuracy and it is independent of the number of approximations.
In contrast to FEM, the suggested method uses essentially fewer equations to solve the structure. Owing to high accuracy of solution of equilibrium equation we need fewer nodes to satisfy boundary conditions. These nodes are placed only at the contour.
There are four main techniques to calculate derivatives in the literature [19]: hand-coded analytical derivatives, numerical differentiation, symbolic differentiation and automatic differentiation. Hand-coded analytical derivatives are exact but time consuming to code, error prone and not automated. Numerical differentiation is easy to code. However, as an estimation approach, it is known to produce inaccurate results. Symbolic differentiation is a technique for computing derivatives of mathematical expressions via symbolic manipulation. This technique is used by computer algebra systems. It is exact, but memory intensive and slow.
In author's computer program automatic differentiation (AD) is used, because it is free from the drawbacks of techniques described above. Unlike symbolic differentiation, which operates on mathematical expressions, automatic differentiation operates on code. Automatic differentiation leverages the chain-rule-based for evaluating the derivatives. It is exact and speed in comparable to hand-coded derivatives. We calculate higher-order derivatives using autograd library. Finally, using jax [20]-the successor of the autograd library-it is possible to run our Python and NumPy programs not only on CPUs (central processing units) but also on much more efficient GPUs (graphics processing units) and TPUs (tensor processing units). Therefore, significant acceleration of calculations can be expected.
The purpose of this paper is construction of effective analytical-numerical method to solve of thin symmetrical plates. The suggested method is based on continuum material model. General solution of equation of equilibrium (1) is presented as sum of general solution of uniform Equation (2) and particular solution of non-uniform Equation (1). They are independent between themselves and contain theoretically infinite number of unknown parameters which allow the performance of the equation of equilibrium and boundary conditions with given accuracy.
General solution of equilibrium equation includes functions of shape of deflection. Functions of shape of displacements, moments, shearing forces and generalized shearing forces have been obtained using automatic differentiation. Such structure of solution allows easily to model various boundary conditions. Obtained relations create the mathematical model of a plate.
Equation of equilibrium is performed by collocation method in each point on the plate surface which are called surface nodes. They are not vertexes of finite element. Improvement of accuracy of solution is achieved by increasing of unknown parameters with author's program. General solution of Equation (2) contains separate unknown coefficients which allow the performance of boundary conditions written in the points at the plate contour. They are called edge nodes. author's program fulfils generation of edge and surface nodes and their distribution.
To obtain solution of the plate with proposed method we at first perform equation of equilibrium in the separate points on the plate surface and obtain particular solution of the problem. Next, we expect that this solution is unchanged and satisfy boundary conditions in the separate points at the plate contour.

Literature Review of the Methods for Solving Plate Bending Problems
Exact solution of the boundary problems of the plate theory are possible only for simplest cases. So far, many methods of approximate solution of plates have been developed. All of them can be divided into two basic groups: analytical [1,2,6,7] and numerical ones [2,9,10,15]. Their analysis is given in fundamental monograph by editor Czeslaw Wozniak [21].
Formulation of the problem of the plate theory usually is performed one of the three techniques: system of equilibrium differential equations, variation approach and method of integral equations.
Solution of differential equations of equilibrium causes development of analytical methods. Numerical methods are less exact but they allow the solving of arbitrary configuration plates [2,10] at the expense of a significant increasing the number of unknowns. It increases the time for calculations and leads to accumulation of computational errors.
Ritz-Timoshenko method is based on the principle of minimum of the structure potential energy (Principle of Virtual Works).
Here Π is potential energy of the structure; V-specific potential energy; A is work of outer load applied to plate structure; δ-Lagrange variation symbol.
Minimalization of error R of approximate solution of Equation (4) with respect to test functions is object of Bubnov-Galerkin method. Navier's method only deals with rectangular plate free supported at the contour and arbitrary loaded. Lévy's method allows the solving of the plates with two opposite edges free supported. Other edges can be arbitrary fixed.
Method of boundary collocation demands satisfaction of integral equations and Trefftz collocation method demands to perform equilibrium equation. In both methods, the boundary conditions are performed at separate points called collocation points at the edges of the plate. Trefftz method does not contain singular integrals and it is more exact in comparison with other methods. Variant of collocation method is worked in paper [36] for a plate inscribed in rectangular contour. Collocation points are Chebyshev nodes. The authors argue that with the number of computational nodes increasing, the accuracy is not improved because of the round-off error effect.
Finite element method (FEM) follows from variation formulation and finite boundary method (FBM) does from boundary integral equation. In the paper [34] four nodal plate element based on Reissner-Mindlin plate theory is worked. Mesh refinement increase the accuracy of the solution. Finite element modeling of stiffened and unstiffened orthotropic plates is presented in the paper [37].
BIM is part of general method BEM associated with calculation singular integrals containing in BEM. According to global element method (GEM) [38] considered area is divided into some simple subareas (greater than in FEM). Division is defined by configuration of the area, solution typologies and mechanical properties of material. In each area approximate functions are taken as linear combinations of polynomial, trigonometric and exponential functions. For example, linear global or superelement has three degrees of freedom. Interpolation functions are taken in the form of second order polynomial. Boundary and continuity conditions at the common edges of the subareas are performed in manner of principle of virtual work. Fault of the method is difficulty of division of considered area into elements and selection of the functional for minimization.
Global-local finite element method (GLFEM) [11] uses Ritz-Timoshenko method at the one part of the plate and finite element method at the another. Fault of this method is necessity of individual treatment of each problem and unstudied effectiveness. According to this method structure is first analyzed as a whole. Therefore, all the local details obtained from the global analysis are considered to be boundary conditions in local details analysis. The global-local analysis is especially useful in the stress concentration analysis, and so on. More information on global-local finite element method can be found in [44][45][46][47][48].
The extended finite element method (G/XFEM) [3,[49][50][51][52][53][54][55] is generalization of classical finite element method and it was especially developed for modeling structural problems with discontinuities. According to that generalized global approximation of the deflection of the plate is presented as [51]: where u j is a nodal parameter associated with finite element shape functions N j (x), b ji are nodal parameters associated with G/XFEM shape functions-N j (x) · L ji (x). L ji are so-called enrichments functions that hold the information about the problem solution.
In the last ten years, there has been a growing interest in meshless methods [56,57]. These methods do not demand of mesh as finite element ones. Area discretization is accomplished by a set of ordered or scattered nodes that are independent of the plate configuration. This approach does not cause difficulties in satisfying boundary conditions. The method is a simply, effective and attractive.
Novel plate and shell analysis methods also include analysis of the Kirchhoff plate using rational Bézier triangles in isogeometric analysis coupled with a feature-preserving automatic meshing algorithm [58] and modeling using a NURBS-based isogeometric analysis (IGA) approach [59].
Some results obtained by analytical and numerical methods are given below. Rectangular uniformly loaded plate free supported at two opposite edges or rigidly clamped at the contour is considered in the paper [60]. To solution the problem Lévy's method was used here. In each concrete case particular solution is taken separately. In similar manner the same plate free supported at the corners was considered in the paper [61].
The exact solution of orthotropic rectangular cantilever arbitrary loaded plate is obtained in [62]. Rectangular plate free supported at the contour was solved [63] using Bubnov-Galerkin variation method. Various kinds of outer load were considered. Using combination of functional series S. Timoshenko obtained original analytical solution of thin orthotropic plate rigidly clamped at the contour [8]. Krjukov using spline approximation method obtained solution of parallelogram and trapeze plates [26].
Trefftz method boundary collocation was used to solve thin anisotropy plate of various geometry [29]. The obtained results are in perfect agreement with known analytical and numerical solutions. Thin rectangular cantilever plate [28] and irregular thin plate on Winkler foundation [36] were also solved with this method.
Finite element method was used [64] to solve the anisotropic elastic plates with holes and cracks. Method of boundary equations is applied to bending problem of the thin plate [65]. Plate is supported on the columns. Substitute forces at the edges and concentrated forces in the corners are not introduced in the given formulation. The problem was solved by boundary integral equations method. Collocation points were placed outside of the plate. This approach allowed to eliminate singular integrals in calculation.
Bending analysis of rectangular moderately thick plates with spline finite element method is subject of the paper [39]. Analysis of deflection of regular and irregular plates with matrix method is fulfilled in the paper [43].
Boundary element method was used to solve stiffened plates [37]. Elastic bending problem of such plates were investigated with boundary element method [66].
Fracture mechanics problems for plane stress and plane strain state solved with extended finite element method is considered in the paper [51]. With help of this method Reissner-Mindlin cracked plate was solved [3]. In article [67] bending by concentrated force of a cantilever strip with a through-thickness crack perpendicular to its axis was solved using linear conjunction method.
The new analytical-numerical method to analysis of symmetric Kirchhoff plate is presented in this paper. The method combines the advantages of an analytical (Navier and Lévy) as well as numerical (Trefftz collocation and finite element) methods. Equation of equilibrium is performed exactly in whole area within the plate. Kinematic and static boundary conditions are satisfied only at the separate nodes at the plate edges with the high accuracy.
The suggested method is close to Trefftz collocation method and meshless ones. It was tested on examples of cantilever plate [68], isotropic irregular plate [69], orthotropic plate resting on Winkler's foundation [70,71], double-connected isotropic rectangular plates [71] and plate connected with truss [72]. Its improved version is given in the presented paper. Polynomial which was introduced in our previous publications for improvement of calculation results has been removed here as unnecessary.

Materials and Methods
Plate is a body constrained by two parallel planes and lateral surface. Upper and bottom planes are bases of the plate. Plane equidistant from the bases is called middle plane. Line of the cross-section middle plane with lateral surface is plate contour. Further we will be designate plate P, surface S and contour C and we take into account the assumptions proposed by Kirchhoff.
Let us consider thin isotropic symmetrical plate with at least two symmetry axes which has thickness h and surface S. Plate is referred to right-handed Cartesian coordinate system Ox 1 x 2 origin at its geometrical center and constrained by contour C (Figure 1a). We assume that it is balanced with outer load q(x 1 , x 2 ) applied to surface S and complied with the boundary conditions. Symmetry of the plate geometry, mechanical properties, boundary conditions and outer load are taken into account in this paper. Plate is defined by coordinates of its vertices We choose the solution of Equation (1) in the form of sum of general solution w 0 of uniform Equation (2) and particular solution w * of Equation (1).
General solution of Equation (2) we take in the following form [73].
Einstein summation rule is used here. According to this convention indices repeated twice in a single term imply the summation is to be done. Summation indices in formula (7) are k, p and s. As particular cases, we obtain Lévy (p = 1, s = 2) and Timoshenko (p = 1, s = 1, 2) forms of representation. In the above p = 1, 2; s = 1, 2; k = 1, . . . , K, where K is several approximations of the solution and determines its accuracy. The K is greater the accuracy of the solution is better.
T kps (x s ) are given trigonometric functions Here Parameters a s (s = 1, 2) are model constants. They must be given a priori. In suggested model they are taken as sizes 2a s (s = 1, 2) of rectangular contour L described on plate contour C (Figure 1b).
Unknown functions f kps (x s ) we take in the form Because indices k, p, s appear at the both sides of this formula, summation according to these indices does not perform.
In the above R kps are unknown coefficients which are determined from boundary conditions at the plate contour, E kps (x s ) are unknown functions. Substituting expressions (7)-(10) into Equation (2) we come to system of two unbounded forth order homogeneous differential equations with respect to these functions [70].
where κ kps = γ ks , p = 1, Functions E kps (x s ) are taken as the following Substituting expression (13) into Equation (11) we come to system of characteristic equations Let λ kpsν , where ν = 1, 2 be their roots. Then functions f kps (x s ) take the form [69] f kps ( In the above Functions B kpsν (x s ) are called the basic functions of the solution.
According to (10) general solution of Equation (2) takes the form Functions are dependent on the solution of the problem and are called shape functions of the plate deflection. General solution of Equation (1) we obtain substituting (17) into expression (6) w(x 1 , Here (1). Functions W * (x 1 , x 2 ) can be present in the form: Functions W mnpq (x 1 , x 2 ) are dependent on outer load and called force functions of the plate deflection. For each specific problem they can be obtain using special procedures.
With the deflection of the plate w and using relations known from theory of thin isotropic plates, we obtain the expressions for tangent displacements [71,72] moments shearing forces and generalized shearing forces In the above functions U kpsν , V kpsν , and so on are partial derivatives of shape functions W kpsν . Similarly functions U * , V * , etc. are partial derivatives of force functions W * . They have been calculated using automatic differentiation [19,74].
Expressions (19)-(24) constitute mathematical model of considered plate. According to this model plate is in internal equilibrium with outer load q(x 1 , x 2 ), but boundary conditions are not satisfied here. Unknown coefficients R kpsν are degrees of freedom of the plate. They are determined satisfying boundary conditions at the separate points of the contour. Let us consider technique of their generation.
Due to the symmetry of the plate, we can consider only its quarter. At first, we introduce two sets X 1 and X 2 of points (x 1m , x 2n ) with numbers (m, n) placed at the positive semi axes of the plate and called them initial points.
Projections of these points onto the contour C are called current nodes and designated K rm and L rn correspondingly, where r is number of the edge (Figure 1b). It is seen that edges parallel to coordinate axes contain current nodes generated by either set X 1 or set X 2 but nodes placed at the inclined edge are generated as set X 1 so X 2 . If two initial points of set X 1 and X 2 fall into one node at the inclined edge we have double current node. Projections of these points into the vertex V i , with number i we call corner nodes and designated K ir , L ir [69]. Union of current and corner nodes we called edge nodes.
Let us consider for simplicity rectangular plate (Figure 2). Its contour has no double nodes. In this case, we can say that each edge node is generated by single initial point. From this reason each function f 1ps (x s ) corresponds to single initial point. At each edge node two boundary conditions are written.
x 1 Corner nodes, say K 32 and L 33 , have the same coordinates (a 1 , a 2 ) but they lie at the different edges of the normals n 2 and n 3 . Thus, they are different nodes. Each edge node is generated by single initial point. At the first approximation (K = 1) we have eight unknown coefficients R 1psν . Thus, eight boundary conditions written at four edge nodes can be satisfied. To that four initial points must be chosen at the semi coordinates axes (x 1m = 1, 2; x 2n = 1, 2). Following for approximation K there must be 2K points x 1m = 1, . . . , 2K; x 2n = 1, . . . , 2K. For illustration of effectiveness of the proposed method as an example we consider rectangular plate ( Figure 2) with thickness h and plane sizes 2a s (s = 1, 2). This example allows easily compares author's results with results obtained by other available methods.
To obtain the particular solution w * of Equation (1) we present outer load q(x 1 , x 2 ) in the form of double trigonometric series [70,73] q(x 1 , For rectangular plates parameters q mnpq can be determined as coefficients of expansion of the outer load in Fourier series on the surface of the plate. New symbols have been introduced here. T m1s (x s ) = cos (γ ms x s ) ; T m2s (x s ) = cos (δ ms x s ) ; For rectangular plates force functions W mnpq are taken in the form Substituting expressions (25), (28) into Equation (1) with taking into account (19,20) and equating coefficients at the same expressions in the both sides of Equation (1) we find coefficients C mnpq expressed by outer load.

Results
Let us consider three kinds of rectangular plates under various boundary conditions:

1.
Rectangular plate with two opposite edges simply supported and the other free, 2.
Rectangular plate simply supported at the contour, 3.
Plate clamped at the contour.
Plates are loaded by particular load where and Here q 0 = 10 kPa is intensity of the load. During the calculations the following geometric and mechanical parameters of model were take into consideration: plate sizes 2a 1 = 8 m, 2a 2 = 4 m, thickness h = 0.2 m. Young's modulus is equal to E = 30 × 10 9 Pa and Poisson ratio ν = 0.2.
The problem has been solved in third approximation (K = 3) using 12 initial points. Consequently, we have 12 edge nodes. Distribution of initial points for (K = 3) approximation is given in table: Due to symmetry of the problem deflection w(x 1 , x 2 ) and the moments M 11 (x 1 , x 2 ), M 22 (x 1 , x 2 ) are symmetrical functions of variable (x 1 ,x 2 ), while slops ϕ 1 (x 1 , x 2 ) and ϕ 2 (x 1 , x 2 ) are antisymmetrical ones.

Rectangular Plate with Two Opposite Edges Simply Supported and the Other Free
Following boundary conditions must be satisfied (see Figure 3) w(x 1 , x 2 )           Bending moment M 11 at the middle cross-section ( Figure 10) is essentially greater than moment M 22 at this section ( Figure 11). Bending moment M 22 and generalized shearing force V 2 are almost zero at the whole simply supported edges (x 2 = ±a 2 ). At the simply supported edges (x 1 = ±a 1 ) deflection and slope ϕ 2 are also zero (values of the order of 1 × 10 −18 and 1 × 10 −17 ). All boundary conditions are performed at the current and corner nodes.
Extreme values of these variables are collected in Table 1. The maximum values are close to the maximum values obtained in the ABAQUS package (Appendix A.1), but the minimum ones for bending moments do not coincide.
Minimum values for the bending moments in the presented examples are values on the edge. Similar observations concerning the comparison of the FEM (ABAQUS) results were included in other works, e.g., in [75,76] where it was indicated that FEM is unable to meet the zero-bending moment at the edge. The explanation of the reason for discrepancy between present method and FEM is included in the Discussion. Compliance with boundary conditions in the midpoints at the edges is presented in Table 2.

Rectangular Plate Simply Supported at the Contour
Following boundary conditions must be satisfied (see Figure 12) w(x 1 , x 2 ) x 1       Extreme values of these magnitudes are collected in Table 3. They exactly coincide with the Timoshenko solution [7]. The maximum values are close to the maximum values obtained in the ABAQUS package (Appendix A.2), but the minimum ones for bending moments do not coincide. Compliance with boundary conditions at the midpoints of the edges is presented in Table 4. In this example numerical studies have been shown that boundary conditions are satisfied with high accuracy − 1 × 10 −30 for deflection and 1 × 10 −20 for bending moments.

Rectangular Plate Clamped at the Contour
Following boundary conditions must be satisfied (see Figure 19) w(x 1 , x 2 )      In this example, numerical studies showed that boundary conditions were met with high accuracy − 1 × 10 −8 for deflection and 1 × 10 −7 for slopes. It is seen that symmetry and antisymmetry conditions are performed. The deflection of the center of considered plate is less than previous plate. At the middle cross-section, the bending moment in the center of the plate is smaller than at the edge.
Results of calculations are presented in Table 5. As in previous examples the maximum values of theoretical results are close to the maximum values of results obtained in the ABAQUS package (Appendix A.3), but the minimum ones for bending moments do not coincide.
The numerical results of bending problem of clamped rectangular thin plates presented in [77] are in excellent agreement between the exact numerical-analytical solution obtained by the generalized integral transform technique (GITT) and FEM obtained by the commercial program ABAQUS, but during the analysis, the plates were discretized using sufficiently refined S4R elements − 200 elements in the x direction and 200 × c elements in the y direction. This means that sufficient accuracy was obtained for 40,000 elements, while in the examples presented in this article, the model consists of 512 elements. Compliance with boundary conditions at the midpoints of the edges is presented in Table 6. Boundary conditions are satisfied exactly. Table 6. Compliance with boundary conditions at the midpoints of the edges.    Plots are constructed for all approximations from 1 to 5. It is seen that beginning from second approximations (K = 2) they are coincided exactly. This confirms the choice of (K = 3) for calculations.

Discussion
Solution of equilibrium equation consists of two parts: general and particular solution which are independent of each other. The general solution identically satisfies differential equation and, in contrast to FEM, satisfies kinematic, static and mixed boundary conditions. The general solution obtained for a rectangular plate is correct for any plate inscribed in a rectangular contour. Particular solution satisfies non-uniform equation of equilibrium using boundary collocation method. It also allows the consideration of plates with arbitrary configurations. Since symmetrical problems are considered, deflections and bending moments are symmetrical functions of variables (x 1 , x 2 ) but slopes are nonsymmetrical ones.
Advantages of suggested method: • Simplicity of structure modeling, • Possibility to define static, kinematic and mixed boundary conditions, • High accuracy and efficiency of calculation, • Meshless approach for solving the problem.
Neither a surface mesh as in FEM nor a boundary mesh as in BEM are required. Introduced nodes are not vertex of finite elements. Surface and edge nodes are generated automatically. Contrary to FEM [78], the presented method is based on continuum model of material. For this reason, operations such as structure discretization and finite element aggregation are unnecessary. Singular plate can be considered to be one global or macro plate element.
Disadvantages of presented method: • It cannot be used directly to solve fracture mechanics problems, • It is worked for simply connected structures; for each additional contour, it is necessary to introduce the new block of basic functions, • It gives good results for regular plates; plates with irregular geometry are solved with less accuracy, • Unlike FEM, system matrix is consisted from zero and non-zero blocks.
Tables 1-6 contain the minimal and maximal values of static and kinematic magnitudes of the plate obtained by author's program based on proposed method written in Python programming language, using finite element method in ABAQUS package, and based on the formulae for a simply supported plate according to Navier's method, included in Timoshenko's monograph [7].
It is seen that analytical and numerical results of kinematic magnitudes coincide. Maximal values of the bending moments are close but they do not coincide. Twisting moments are differ significantly. Numerical values of bending moments are not zero along the edges of the plate. It follows that static boundary conditions are not met in ABAQUS package.
The reason for this discrepancy is as follows: • Presented method is based on Kirchhoff theory of the plates which is deformation theory of zeroth-order where deformations of a cross-section are omitted. Instead generalized shearing force as a sum of shearing forces and derivative of twisting moments is introduced. • ABAQUS package is based on Mindlin plate theory which is deformation theory of first-order where deformations of a cross-section are taken into account. Twisting moments and shearing forces are independent ones. It is reason that obtained values of bending and twisting moments are different.

•
Relatively small number of finite elements. The Student Edition of the ABAQUS package is restricted to 1000 nodes.
In analytical approach kinematic and static boundary conditions are performed directly at the separate points at the plate contour. Numerical methods perfectly satisfy kinematic continuity equations and kinematic boundary conditions only at the separate nodes. Equations of equilibrium and static boundary conditions are fulfilled approximately in the whole area of the plate using Principle of Virtual Work. It can be reason that kinematic boundary conditions are satisfied very well but static ones are performed with less accuracy in Finite Element Method.

Conclusions
There are many analytical, numerical and analytical-numerical methods to solve thin plates. Most of them are based on the theory of thin plates (see Section 2). Method presented in the paper is one of them. The essence of the method is as follows: • to obtain with high accuracy particular solution of equilibrium equations at the separate surface nodes using boundary collocation method, • to generate set of the initial points at the coordinate axis and then their distribution at the plate contour, • to write boundary conditions at each edge nodes; number of these nodes always corresponds to number of unknown parameter of the model, • to solve system of boundary equations, • to calculate the desired magnitudes, • to prepare the results.
The major findings of this research are summarized in the following: • It is shown that chosen direction is actual and important from a practical point of view. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following notations are used in this manuscript: Twisting moment per unit length of sections of a plate perpendicular to x 1 axis Q 1 , Q 2 Shearing forces parallel to x 3 axis per unit length of sections of a plate perpendicular to x 1 and x 2 axes, respectively V 1 , V 2 Generalized shearing forces parallel to x 3 axis per unit length of sections of a plate perpendicular to x 1 and x 2 axes Figure A1. The deflection of a rectangular plate with two opposite edges simply supported and the other free solved using FEM.  Figure A2. The bending moment M 11 of a rectangular plate with two opposite edges simply supported and the other free solved using FEM.  Figure A4. The deflection of a simply supported rectangular plate solved using FEM.  Figure A5. The bending moment M 11 of a simply supported rectangular plate solved using FEM.  Figure A7. The deflection of a rectangular plate with clamped edges solved using FEM.  Figure A8. The bending moment M 11 of a rectangular plate with clamped edges solved using FEM. return ((x(s)(x_1, x_2) / a(s)) * (np.sinh(kappa(k, p, 3-s) * x(s)(x_1, x_2)))) return wrapper