Copyright © 2022 O.V. Bagrii. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Problems of Tribology, V. 27, No 2/104-2022, 104-111 Problems of Tribology Website: http://tribology.khnu.km.ua/index.php/ProbTrib E-mail: tribosenator@gmail.com DOI: https://doi.org/10.31891/2079-1372-2022-104-2-104-111 Plane problem of discrete environment mechanics O.V. Bagrii Khmelnytskyi National University, Ukraine E-mail: avadaro@yahoo.com Received: 28 April 2022: Revised: 28 May 2022: Accept: 25 June 2022 Abstract Many engineering problems related to the design of structures and machines, the mathematical description of technological processes, etc., are reduced to the need to solve a plane problem for materials with a significant effect of internal friction on their deformation. Such materials include a large class of materials in which the compressive strength is greater than tensile. These are composite materials, concretes, rocks, soils, granular, loose, highly fractured materials, as well as structurally heterogeneous materials in which rigid and strong particles are interconnected by weaker layers. The laws of deformation and destruction of such materials differ significantly from elastic ones. A feature of these laws is an increase in resistance to shear deformations and an increase in the strength of materials with an increase in the magnitude of compressive stresses. This is associated with the influence of internal Coulomb friction on the process of their deformation in the limiting and boundary stages. The need to formulate and solve a special boundary value problem for materials with significant internal friction is because the results of solving problems using models of elasticity and plasticity differ significantly from experimental data. The difference increases when approaching the limiting state of discrete materials and depends significantly on the structure of the material and operating conditions. The boundary value problem of the mechanics of a deformable solid is formulated as a system of equations of three types: static, geometric, and physical. For all linear and physically nonlinear problems, provided the deformations are small, the first two groups of equations remain the same. Thus, these differences can be attributed to the inconsistency of the accepted in the calculations of physical relations "stress - strain" and the real laws of deformation of these materials, which are more complex rheological objects than structurally homogeneous solids, liquids or gases. The article uses an approach where the material is immediately considered as quasi-continuous, and the physical equations are based on the experimentally obtained relationships between the invariants of the stress and strain tensors, which consider the influence of both molecular connectivity and internal Coulomb friction. Key words: Plane problem, internal friction, compressive stress, variable deformation modules, nonlinear physical equations Introduction Plane problem for materials with significant internal Coulomb friction is formulated as a boundary value problem for a flat, inhomogeneous, physically nonlinear area filled with a material, the deformations of the forms of change of which are significantly affected by the values of compressive stresses [1]. The problem is to determine the stress and strain fields when the region is perturbed by force or kinematic factors The material in the calculation area is considered quasi-continuous. The material deformation model is a combination of a perfectly cohesive material model (Prandtl model) and an internal friction model (Coulomb model). The physical relations of the model are written in the form of relations of mechanics of deformable solids but with variable modules of deformation, the values of which depend on the achieved level of stresses and http://creativecommons.org/licenses/by/3.0/ http://creativecommons.org/licenses/by/3.0/ http://tribology.khnu.km.ua/index.php/ProbTrib https://doi.org/10.31891/2079-1372-2022-104-2-94-103 Problems of Tribology 105 strains and are determined by testing macro samples of material under flat deformation. A feature of the formulated problem is using physical equations for the active deformation process with variable deformation parameters that are different at each point of the computational area. Therefore, the problem's solution can be obtained only by specially developed iterative procedures. The analysis of the known numerical methods allowed to choose finite elements with some limitations as the basic method for the realization of the formulated physically nonlinear problem of a flat inhomogeneous area. First, the fulfillment of the condition of invariance of stresses within one element, which is achieved by a special choice of the shape of the finite element and approximating functions. This makes it possible to assign values of variable deformation modules at a particular stage of calculation not for each point of the calculation area, but for each finite element, which allows the implementation of special iterative algorithms for solving the problem by numerous methods. The second fundamental condition is the need to consider the volume forces of gravity. Without this, the discrete material may not perceive an external load. Therefore, it is necessary to formulate a nonlinear mathematical model of the environment using physical relations that consider the influence of internal Coulomb friction on the deformation of the environment under plane deformation conditions. Mathematical formulation of a plane problem To formulate the problem, we take from the array of the medium operating under conditions of plane deformation, a flat design region О of unit thickness 1h (Fig. 1), which can have finite dimensions or be considered as part of an infinite plane. On some section S(R) of the contour, external loads (force boundary conditions) can be specified, and known displacements (kinematic boundary conditions) on the other S(u). On the part S(0), the boundaries of the displacement area may be zero. The laws of material deformation are described by nonlinear dependencies that consider the influence of internal friction. In this formulation, the problem is reduced to determining the fields of stresses, strains and displacements (σ(x, y), ε(x, y), u(x, y)) of the computational area, which correspond to the accepted model of the environment and given boundary conditions. The problem is formulated by systems of equations: equilibrium equations: ;                    y xyy x xyx V xy V yx (1) geometric Cauchy relations: x ux x    ; y u y y    ; y u x u xy xy       ; x u y u yx yx       ; (2) nonlinear physical equations of discrete materials: Fig.1. Calculation area of a flat problem 106 Problems of Tribology .; 2 2 2 2 ; 2 2 2 2                    yхзмxy xyзмxy y змзм x змзм y y змзм x змзм x γGτ γGτ ε GK ε GK σ ε GK ε GK σ (3) The last relations, as noted earlier, differ from the equations of the plane problem of the theory of elasticity in that the modules змK and змG are functions of the stress-strain state and not constants of the material.. Mathematical formulation of the plane problem of the mechanics of discrete materials in a convenient for numerical implementation of the matrix form can be represented as follows. For a flat region О, on the boundary of which force or kinematic boundary conditions are specified, determine the stresses     yx, , deformation     yx, and displacement     yxuu , corresponding to the calculation scheme and the accepted specific laws of deformation of the region's material. The problem is reduced to solving a system of matrix equations:      VB T  – differential equilibrium equations (4);     uB – differential geometric equations (5);      D – physical equations for discrete material (6) considering the experimentally obtained invariant nonlinear physical relations to determine the variable parameters змG , змK : ,P m n Gзм      1 1 2 змзм GK , (7) and boundary conditions:      S T RA  – in S(R);    Suu  – in S(u); (8)    0u – in S(0), where  SR – vector of known forces acting on the boundary S(R),  TA – direction cosine matrix   ; 0 0        ml lm A T (9)  Su – known displacements at the boundary S(u);  0 – lack of movement at the boundary S(0). In addition to the formulated conditions for discrete materials (primarily for bulk, granular and grainy) a fundamental restriction is introduced on the impossibility of tensile stresses. If we introduce the rule of signs, according to which the compressive stresses are considered positive, this restriction is formulated as .02  (10) A feature of the formulated plane problem is that the physical equations do not include the elastic characteristics of the material G and K. Instead, functions are used with different values at each point of the computational area, which depend on the level of stress-strain state reached in it. Therefore, to solve the formulated problem, specially developed numerous iterative procedures are used. There are many effective numerous methods that could become the basis for the development of iterative algorithms necessary to solve the problem. The most famous of them are: the finite difference method (FDM), the finite element method (FEM), the boundary element method (BEM). The analysis of the features and capabilities of these methods made it possible to choose the implementation of the formulated physically nonlinear problem of a flat inhomogeneous region as the basic finite element method with some restrictions. Finite element formulation of a plane problem The idea of the method is to replace the computational area with its discrete model formed by a system of finite elements that interact with each other in nodes. At these points, the continuity conditions are always Problems of Tribology 107 provided. Each element is deformed according to the laws of deformation of the material of the calculation area. The distribution of stresses and strains inside the element depends on its shape and the form of the accepted approximating functions. The use of the FEM for continuum areas reduces the solution of the boundary value problem to the problem of structural mechanics, which is solved by the known method of displacements. The convenience of using FEM for solving nonlinear problems is due to the simplicity and clarity of interpretation of all stages of the calculation, the possibility of step-by-step control of the results. However, the most important feature of the method is that for specially selected elements of the element and the approximating functions of stress and strain will not change within one element. Since the modulus of deformation of discrete materials depends on the achieved level of stress, the values of the modulus of strain at a particular stage of the calculation in this case can be assigned not for each point of the calculation area, but for each finite element. This makes it possible to implement special iterative algorithms for calculations by numerical methods. The sequence of solving a plane problem based on the finite element method is quite fully described in the scientific and methodological literature [2, 3] and is reduced to the following operations: 1. Create a discrete model of the computational area - divide the solid domain into a grid of finite elements (see Fig. 2). 2. External loads at the boundary of the region S(R), and body forces lead to forces applied at the nodes  R , and known displacements on the boundary S(u) lead to displacements of the nodes   . The condition of the impossibility of displacements on the boundary S(0) is satisfied by introducing appropriate restraints at the nodes on the boundary of the region. 3. For each finite element, a stiffness matrix [ke] is formulated that connects the vectors of nodal forces and displacements of nodes belonging to the same element. 4. Formulate a global stiffness matrix [K] of the system of elements, consider the "contribution" of each finite element to the stiffness of the entire system. 5. Formulate and solve a system of canonical equations of structural mechanics, most often the displacement method, as a result, nodal displacements are found, and through them - strains and stresses in finite elements. { { Fig. 2. Discrete model of computational area Most of the described steps are standard for FEM. However, to implement the previously formulated physically nonlinear problem, it is necessary to consider the following features. 1. The shape of the finite element and the adopted approximating functions must ensure at each stage of calculation the invariability of the achieved values of stresses and strains within one element, which allows assigning variables of deformation modules not for each point of the region but for each finite element. This is achieved by choosing a basic finite simplex element in the form of a triangular prism of constant length L = 1 = const (Fig. 3). Fig. 3. Simplex element for a plane problem 108 Problems of Tribology An approximating function that describes the nature of changes in displacements from the coordinates of x, y points is assumed to be linear. The displacements {u} of an arbitrary point in the middle of an element are expressed in terms of the nodal displacements {δe} of its nodes by the following relation:        ,emjie y x INININN u u u           (11) where I – single dimensional matrix 2×2; N – shape function , 2 A cybxa N   (12) where a, b, c - coefficients determined by nodal displacements; A – element area (volume 1 AV ). The deformations εx, εy, γxy of the points of a finite element according to the Cauchy relations are derived derivatives of the displacements {u}. For the accepted linear approximation, the derivatives do not depend on the x, y coordinates of the points of the element. Therefore, the deformations in the middle of each element are constant. This can be easily verified by performing the following differential operations.        ,,,                             m j i mjiee xy y x e BBBB (13) where mji BBB ,, – submatrix of the matrix  eB of the differential operator. For example, the submatrix  iB looks like:   .0 0 2 1 0 0                                       ii i i i bc c b A x N y N y N x N B (14) Similarly written matrices  jB ,  mB . Thus, the matrix  eB , and therefore the deformations  e do not depend on the coordinates of the points of the element There is an unambiguous relationship between deformations  e and stresses  e :      ,eзмee D  (15) where   змe D – matrix of deformation parameters. Therefore, the values of the components of the stress tensor within the finite element will also be constant at each stage of the load. Therefore, the choice of a finite simplex element in the form of a triangular prism with a linear displacement function satisfies the above requirements. 2. The second fundamental feature of the problem under consideration is the need to take into account the volume forces of gravity. Without this, the computational area of a discrete material sometimes cannot perceive an external load. As shown by O. Zenkevich [4], body forces X, Y can be reduced to nodal forces using the integral relation     .        dxdy Y X NR T e (16) For the introduced shapes of the elements and the linear approximating function, after integration, we obtain that the force of gravity is equally distributed between the three nodes of the element. The procedure for reducing the known values of the distributed load P and displacements u at the Problems of Tribology 109 boundary of the computational domain to nodal forces {R} and displacements {δ} does not differ from that adopted for linear problems. 3. The next stage of calculation - the formulation of the stiffness matrix of the element - under the accepted assumptions is also simplified. In the general case, the stiffness matrix of a finite element for a plane problem is determined by the expression:        . dxdyBDBk ee T ee (17) Considering the above assumptions, expression (4.22) takes the form:         .ABDBk eзмe T ee  (18) where A – the cross-sectional area of the finite element plane xy . The components of the matrix  eB within the finite element do not depend on the coordinates, and the matrix of deformation parameters   змe D contains variable modules Gзм, Kзм, the values of which depend on the stresses and strains achieved in the center of gravity of the finite element.   , 200 022 022 2 1              зм змзмзмзм змзмзмзм змe G GKGK GKGK D (19) where ,P m n Gзм   . 1 1 2    змзм GK The formation of the global stiffness matrix [K] of the system of elements of the discrete model of the computational area is a well-known procedure for “assembling” stiffness matrices of elements, described by O. Zenkevich [4]. After the formation of the global stiffness matrix [K], the continuous calculation area is considered as a normal mechanical system of deformable elements, the calculation of which is a known method of displacement of structural mechanics and is reduced to finding displacements and forces of nodes. To do this, additional restraints are set in the hinges that connect the finite elements, which make the system kinematically definable. These restraints allow many possible movements of nodes. The principle of minimizing the potential energy of the system (Lagrange principle) is used to obtain the actual displacements of the nodes [5]. The condition of the minimum of potential energy in the method of displacements is written in the form of a system of canonical equations:     ,RK  (20) where  R – vector of external forces applied in the nodes of the system;   – vector of movements of nodes. In their physical essence, canonical equations (20) are equilibrium equations of nodes in which node displacements   are unknown. To obtain the solution of a specific problem, equations (20) must satisfy the boundary conditions. The form of recording these conditions depends on the structure of the basic equations and the known values of forces and displacements at the boundary of the calculation area. If force boundary conditions in the form of known nodal forces are set at the boundary of the calculation domain, they are automatically considered by the vector  R of external forces (the corresponding components of the vector are equated to known forces). For nodes at the load-free boundary, the vector components  R are zero. If the boundary conditions are given in the form of known node displacements at the boundary of the region, the vectors   and  R , as well as the stiffness matrix [K] are partially modified by a known procedure, for example, the procedure described in [4]. The final finite element formulation of the problem is reduced to the following operations. Create a discrete model of the calculation area. To do this, the solid region is replaced by a grid of finite elements in the form of triangular prisms connected in nodes. The necessary matrices of connection of numbers of elements and coordinates of knots, vectors of known loadings and movements, force and kinematic boundary conditions are formed. For each finite element, a stiffness matrix [ke] is formed, which satisfies the conditions of equilibrium and 110 Problems of Tribology continuity, as well as the physical relationship between stresses and strains. These conditions are described by the following equations:      VB T  – equilibrium equation (4);     uB – Cauchy equation (5);       змD – physical nonlinear equations of discrete materials (6). Finite element stiffness matrix         .ABDBk eзмe T ee  (21) According to the principle of "assembly" of the stiffness matrices of the elements, the global stiffness matrix [K] of the computational area is formed and the system of linear equations (20) of the displacement method is written. Solving the system considering the boundary conditions, we obtain a vector of nodal displacements:      ,1 RK  (22) and through it deformations are calculated      B and stress       змD in the middle of the element. The described procedure can be implemented for a physically nonlinear problem only after the development of special computational iterative algorithms. Its essence is that a complex physically nonlinear problem is solved in stages. At each stage, a linear problem is solved in which the values of the deformation parameters of the element are assigned depending on the level of stresses and strains achieved in it at the previous stage. Conclusions In the mathematical formulation of the problem, in addition to the known equations of equilibrium and continuity of deformations, specific physical equations are used, in which instead of elastic constants variable parameters are introduced that depend on the values of stresses and strains achieved at each stage. The described plane physically nonlinear problem of inhomogeneous region is solved by means of specially developed iterative procedures. The finite element method is chosen as the base in these procedures. The choice of finite elements in the form of triangular prisms with a linear function of the form is substantiated. The formulated finite-element matrix relations allow to obtain the solution of the problem by the numerical method using special iterative procedures. References 1. Kovtun VV Formulation of the problem of contact interaction of structural elements with a discrete environment / VV Kovtun, OV Kolesnikova // Problems of modern engineering technology: 4th interuniversity. scientific-theoretical Conf., January 14–15, 2003: Coll. Science. work. - I., 2003. - Part II (special issue). - № 24. - P. 127–130. [in Russian]. 2. Bezukhov N. I., Luzhin O. V. Using methods of the theory of elasticity and plasticity for solving engineering problems. - M. : Higher school, 1974. - 200 p. [in Russian]. 3. Handbook of the theory of elasticity (for civil engineers) / [ed. P. M. Varvak, A. F. Ryabov]. - K. : Budivelnik, 1971. - 418 p. [in Russian]. 4. O. C. Zienkiewicz. The Finite Element Method In Engineering Science / O. C. Zienkiewicz – McGraw- Hill; 1st Edition (January 1, 1971). – 521 p. 5. Bugrov A. K., Grebnev K. K. Numerical solution of physically nonlinear problems for soil foundations // Foundations, foundations and soil mechanics. - 1977. - No. 3. - P. 39-42. [in Russian]. Problems of Tribology 111 Багрій О.В. Плоска задача механіки дискретного середовища Багато інженерних задач, що пов’язані з проектуванням споруд і машин, математичним описанням технологічних процесів та ін., зводяться до необхідності розв’язання плоскої задачі для матеріалів з суттєвим впливом внутрішнього тертя на їх деформування. До таких матеріалів відносять великий клас матеріалів, у яких міцність на стиск більша ніж на розтяг. Закони деформування та руйнування таких матеріалів суттєво відрізняються від пружних. Особливістю цих законів є збільшення опору деформаціям зсуву і збільшення міцності матеріалів з ростом величини стискаючих напружень. Це асоціюється з впливом внутрішнього кулонового тертя на процес їх деформування в дограничній і граничній стадіях. Крайова задача механіки деформівного твердого тіла формулюється як система рівнянь трьох видів: статичних, геометричних і фізичних. Оскільки для усіх лінійних та фізично нелінійних задач за умови малості деформацій перші дві групи рівнянь залишаються тими ж, відмічені розбіжності можна віднести за рахунок невідповідності прийнятих в розрахунках фізичних співвідношень „напруження – деформації” і реальних законів деформування вказаних матеріалів, які є більш складними реологічними об’єктами ніж структурно однорідні тверді тіла, рідини чи гази. В статті використано підхід, коли матеріал одразу розглядається як квазісуцільний, а фізичні рівняння ґрунтуються на одержаних дослідним шляхом співвідношеннях між інваріантами тензорів напружень і деформацій, які враховують вплив як молекулярної зв’язності, так і внутрішнього кулонового тертя. Ключові слова: плоска задача, внутрішнє тертя, стискуюче напруження, змінні модулі деформації, нелінійні фізичні рівняння