Nonlinear analysis of a crossing- beams system on an elastic with usage of the Mathematica software

In this paper, the authors consider the method of calculating an infinite system of cross beams on an elastic base by the variationaldifference method. The system of cross beams on an elastic base is most often modeled as shallow strip foundations for buildings of various functional purposes. The variation-difference method is one of the numerical and analytical methods for calculating building structures, it is based on the variational principles of the Ritz-Timoshenko method and on the minimum of the total potential energy of the entire system according to the Lagrange principle, and is also close to the real operating conditions of the foundation – base. A single-layer isotropic artificial base was used as an elastic base in the work, as an elastic layer limited in thickness. The algorithm of nonlinear calculation is based on the use of the iterative method of elastic solutions. The physical nonlinearity of the material of reinforced concrete beams is taken into account through the asymptotic dependence "Moment-curvature". Numerical approbation of the results of elastic and nonlinear calculations of the system of cross beams on an elastic base was carried out using the MATHEMATICA software package.


Introduction
At the first stages of the development of the theory of beams on an elastic half-space, it was assumed that there should be no special difference in the calculation results in the conditions of spatial and plane problems (like the calculation according to the Winkler hypothesis) [1]. However, further studies [2] showed that the difference in the results is very large and when calculating tape foundations in the conditions of a flat problem, completely unacceptable strength reserves are obtained (see the fourth chapter, Fig. 115, 126, 131 [1]). Problems about beams in spatial conditions are more complicated, and therefore much fewer methods for solving them have been proposed.
However, even Wieghardt in his work [3], based on the Boussinesq formula, pointed out that if the desired law of reactive pressures p (x, y) depends on the coordinates of the element of these pressures, then it is associated with soil precipitation under the beam W (x, y) (formula 32 [1]). And to simplify it, Wieghardt proposed to assume that the reactive pressures are distributed evenly in the transverse direction and, if the soil sediment and the deflections of the beam are equal, take the transversely averaged precipitation.
Proctor [4], as well as Wieghardt, proceeded from the equations (32), (33) [1] and the assumption that the reactive pressures are evenly distributed in the transverse direction, and as the value of the sediment in the equation of equality of the deflections of the beam and the soil precipitation (33) took their value along the middle axis of the beam. As for the deflections of the beam, the beam is so rigid in the transverse direction that its bending in this direction can be neglected. However, apart from assumptions and hypotheses, neither Wieghardt nor Proctor were able to provide a final solution to the problem of calculating an infinite beam in spatial conditions.
Research on solving this problem was continued by V. I. Kuznetsov, who believed that Proctor's solution is rather cumbersome, and the accuracy is low, and in his work [5] proposed another version of Proctor's solution, taking reactive pressures as an unknown function, and not precipitation, as in Proctor's. Kuznetsov's solution is also very timeconsuming and it has not received significant distribution.
The method of B. N. Jemochkin [6], similar to his method for calculating bands in the conditions of a flat problem, turned out to be simpler, but no less accurate. The only difference is that the Boussinesq formula for the spatial problem is used here instead of the Flemann formula for the plane problem. The beam is divided longitudinally into a series of rectangles, inside which the pressures transmitted by the sole of the beam to the ground are considered constant. Stress intensities in each step of the diagram are considered to be unknown. Just like Proctor, Jemochkin neglected the fact that with a uniform distribution of pressures in the transverse direction, a "hole" is obtained and therefore the deflections of the beam and the soil precipitation under it along the longitudinal axis and outside it will not coincide.
This inaccuracy is absent in the solutions of problems about beams of finite length on an elastic half-space [1,2], in which M. I. Gorbunov-Posadov used the same method that he used to solve the strips in the plane problem [7,8]. However, in the case of solving the spatial problem of the contact interaction of infinite strips, ribbons, beams and an elastic base, Gorbunov-Posadov took a double power series with unknown coefficients as the initial equation for the reactive pressures p (x, y), as in the expression for the soil precipitation W (x, y). Moreover, the unknown power series coefficients for soil precipitation A 2i,2k are linear functions of the power series coefficients for reactive pressures а 2i,2k .
The solution of contact interaction problems for flexible structures on an elastic base by the methods of elasticity theory [9] and structural mechanics [10] was further developed in the works of S. V. Bosakov, S. D. Semenyuk, O. V. Kozunova [11][12][13][14][15][16], which took into account the inhomogeneity (layering) of the elastic base, its physical nonlinearity, creep of concrete, etc. complicating parameters of the contacting bodies.
Due to the complexity of solving contact problems, especially for bendable structures, the scientific literature on the use of variational methods in solving contact problems of the theory of elasticity is very scarce. When calculating the system of cross beams on an elastic base, the authors assume that the system of cross beams is a set of orthogonal rigidly interconnected rods located on an elastic base with axes located in the same plane coinciding with one of the main axes of inertia.
The physical nonlinearity of the material of reinforced concrete beams is taken into account through the asymptotic dependence "Moment-curvature". The "Moment-curvature" dependence in the nonlinear calculation of reinforced concrete beams operating under flat bending conditions comprehensively takes into account the nonlinear properties of concrete and reinforcement, anisotropy and inhomogeneity, cracking of the beam material.
It should be noted that for the first time Murashov V. N. proposed to use this dependence in the theory of reinforced concrete [17]. In the works of V. I. Solomin and his students [18], this dependence was successfully applied to the calculation of reinforced concrete structures on an elastic base.
Hypotheses and assumptions. When calculating a system of cross beams on an elastic base under the conditions of a spatial problem, the following hypotheses and assumptions are introduced: • for the calculation area of an elastic base, hypotheses and assumptions of the linear theory of elasticity and small elastic -plastic deformations of Ilyushin are introduced, taking into account the physical nonlinearity of the beam material [9,10,19]; • for the contact zone of the "foundation -base" system -both compressive and tensile stresses can occur between the infinite system of cross beams and the base and there are no friction forces; • the distribution of normal reactive pressures over the width of each beam is considered constant [1]; • the hypotheses and assumptions of the flat bending of the beam are valid for the beam [10].

Problem statement
An infinite system of cross beams on an elastic base of constant bending stiffness EJ х , EJ y under the action of a symmetric load is considered. Due to the symmetry, an infinite system of cross beams is divided into a number of basic fragments in the form of cross-shaped intersections of these beams (Figure 1), connected to the system. And the infinite system is replaced by a set of intersecting beams freely resting on an elastic base.
The linear dimensions of the beams are indicated as l x , l y . The cross-sections of the beams are assumed to be constant. The external load acts perpendicular to the plane of the axes of the system of cross beams. Boundary conditions of the problem. At the boundaries of the accepted computational domain, the horizontal displacements u = 0, ʋ = 0. In the contact zone, the equality of the base sediment to the deflections of the beams is valid.

The algorithm of linear calculation
In a linear (elastic) calculation, the elastic base is replaced by a computational domain for solving a spatial problem ( Figure 2). The base is approximated by a symmetric three-dimensional centre grid with constant steps along the axes: Δх, Δу, Δz. As a result, 96 cells and 175 nodal points were obtained.  The solution of the problem is constructed in displacements, taking the vectors of nodal displacements u i (x,y,z), v i (x,y,z), w i (x,y,z) as unknown components. When solving this problem, the strain energy is calculated for each cell, and then summed up by the volume of the elastic base. In this case, the system of differential equations after replacing the integrodifferential expressions of the energy functionals with finite-difference approximations is transformed into a system of linear algebraic equations, the solution of which allows us to find unknown components of the displacement vector u i (x,y,z), v i (x,y,z), w i (x,y,z).
In order to find the strain energy on the centre of the volume cell of the computational domain, it is first necessary to find the functional of the strain energy of the elastic base for the centres of the cell faces through the known dependencies of the plane problem of elasticity theory: the Cauchy relations and the generalized Hooke law. Since the cell of the computational domain is a parallelepiped with pairwise equal faces, it is sufficient to determine the strain energy for three faces.
We consider the first face of the cubic cell of the finite difference method with the centre at the point k 1 (Figure 4). We write the expressions for deformations ; We consider the second face of the cubic cell of the finite difference method with the centre at the point k 2 ( Figure 5).  ; We consider the third face of the cubic cell of the finite difference method with the centre at the point k 3 (Figure 6). We write the expressions for deformations  ; The energy of deformations of a rectangular cell with dimensions (хy) centered at the point k 1 , averaged will be equal to (4) where , kk E  − the elastic constants on the centre of the cubic cell of the base; i, j − the number of the nodal point along the X and Y axes, respectively.
Similarly to formula (4), the deformation energy of a rectangular cell with dimensions (хz) centered at the point k 2 will be averaged equal to (5) where i, t − the number of the nodal point along the X and Z axes, respectively. The energy of deformations of a rectangular cell with dimensions (yz) centered at the point k 3 is similar to the formulas (4,5), averaged will be equal to (6) where j, t − the number of the nodal point along the Y and Z axes, respectively The sequence of calculation stages is given by the algorithm of linear calculation of the system of cross beams by the Ritz-Timoshenko method [15,19].
In the works of one of the authors [16], the proposed form of the functional of the total potential energy was previously given, taking into account the deformation energy of the elastic base, which is physically nonlinear and inhomogeneous.
The functional of the deformation energy of an elastic base in a unit volume [19] can be represented as: where E, µ  the elastic constants of the elastic base. Denoting the volume element by dv, the total deformation energy of the elastic base will have the form When composing the functional of the deformation energy of the elastic base, the work of the forces of the elastic base's own weight is not taken into account, since the forces of the elastic base's own weight are balanced by the initial stress state already in the elastic base, and the work of the self-balanced system of forces at small possible displacements is zero. This means that when searching for the complete stress state of the problem under consideration, it is necessary to impose a stress state on the resulting solution from the forces of the base's own weight.
The bending energy of two cross beams is determined by the formula where EJ x , EJ y  the bending stiffness of beams. The deformation energy of the structure is usually identified with the bending energy of the structure, ignoring the shear deformations [10,15]. This is quite justified for the considered infinite system of cross beams.
The external load potential is determined from the following formula The functional of the total energy has the Э U П     (11) Unknown displacements u i (x,y,z), v i (x,y,z), w i (x,y,z) can be found from the condition of vanishing the derivatives of the total energy for each of the displacements, since in a state of static equilibrium, the functional of the total energy Э must have a minimum, that is 0, 0, 0, 1, 2,3,..., , (12) where N − the number of nodal points of the parallelepiped. Numerical approbation of the calculation results for an elastic base is carried out using the MATHEMATICA computer algebra software package [20].
Taking into account the physical nonlinearity of the beam material. In the proposed work, the "Moment-curvature" dependence is approximated in the form of a hyperbolic tangent (13), conditionally shown in Figure 7. We will consider the beam material to be non-linearly elastic where lim M  the maximum bending moment perceived by the beam cross-section; 0 B  initial beam stiffness; 1 ρ  the curvature in this section of the beam. This dependence was widely used by one of the authors in his works to take into account the physical nonlinearity of the elastic medium when studying the stress-strain state of inhomogeneous foundations of foundation structures [21], as well as articulated beams on an elastic base [22]. In the nonlinear calculation, we will use the secant bending stiffness of the beams (see Figure 7).

The results of the linear calculation
To find the solution numerically, the following initial data are taken: l x =l y =4m, l z =6m, x=y=z=1m, EJ х =2000 kNm 2 , EJ y =2000 kNm 2 , E=3,0610 10 MPa; µ=1/6, P=2000 kN.  The internal forces in the beams can be determined by the following differential dependencies, using finite differences

Conclusion
In this paper, the authors propose a variational-difference method to investigate the parameters of the stress-strain state of tape foundations on an elastic base, as an infinite system of cross beams located on an elastic half-space with a limited depth of the compressible thickness. An elastic calculation algorithm is constructed taking into account the physical nonlinearity of the beam material, a program is compiled using the MATHEMATICA computer package, its approbation is carried out and the way to solve the nonlinear problem is indicated. As a result of the conducted studies, it was noticed that under the boundary conditions and numerical data taken in the task, the beams bend in a wave-like manner, which is assumed in full-scale conditions for the selected dimensions of the calculated area and requires clarification of the depth of the compressible thickness of the elastic base.