NUMERICAL METHODS FOR CONTACT ANALYSIS OF COMPLEX-SHAPED BODIES WITH ACCOUNT FOR NON-LINEAR INTERFACE LAYERS

A large number of engineering constructions contain elements that are in the conditions of contact force and kinematic conjugation. For instance, it can be the modified working surfaces of gears, roller bearings, etc. Traditional methods for modeling of contact interaction lead to significant errors in the obtained results or to overly cumbersome numerical models. Another important factor is the lack of adequate modeling of the contact conditions at the boundaries of the bodies. For instance, the impenetration condition for smooth bodies in linearized form is most commonly 10. Ivanyk I. H., Vikhot S. I., Vybranets Yu. Yu. Metodyka rehuliuvannia zusyl v kombinovanykh statychno nevyznachenykh stalezalizobetonnykh konstruktsiyakh // Mekhanika i fizyka ruinuvannia budivelnykh materialiv ta konstruktsiy. 2007. Issue 7. P. 443–453. 11. Teoretychni doslidzhennia napruzheno-deformovanoho stanu kombinovanykh statychno nevyznachenykh metalevykh konstruktsiy / Ivanyk I. H., Vikhot S. I., Pozhar R. S., Vybranets Yu. Yu. // Visnyk natsionalnoho universytetu Lvivska politekhnika «Teoriya i praktyka budivnytstva». 2008. Issue 627. P. 106–111. 12. Prostorovyi rozrakhunok kombinovanykh stalezalizobetonnykh system / Ivanyk I. H., Vikhot S. I., Vybranets Yu. Yu., Ivanyk Yu. I. // Haluzeve mashynobuduvannia, budivnytstvo. 2014. Issue 3 (42). P. 86–91. 13. Eksperymentalni doslidzhennia deformovanoho stanu kombinovanykh statychno nevyznachenykh stalezalizobetonnykh konstruktsiy / Ivanyk I. H., Vikhot S. I., Pozhar R. S., Vybranets Yu. Yu. // Visnyk Natsionalnoho universytetu Lvivska politekhnika. 2007. Issue 600. P. 142–147. 14. Wight J. K., MacGregor J. G. Reinforced Concrete: Mechanics аnd Design. 6th ed. Wilsons Travels Stock, 2011. 1177 p. 15. ACI Innovation Task Group 4. ITG-4.3R-07 Report on Structural Design & Detailing for High Strength Concrete in Moderate to High Seismic Applications. ACI Committee, 2007. 62 p. 16. Kaar P. H., Capell H. T., Hanson N. W. Stress-Strain Characteristics of High-Strength Concrete. Douglas McHenry International Symposium on ConcreteStructures, ACI Publication SP-55, American Concrete Institute. Detroit: MI, 1978. P. 161–185.


Introduction
A large number of engineering constructions contain elements that are in the conditions of contact force and kinematic conjugation.For instance, it can be the modified working surfaces of gears, roller bearings, etc.
Traditional methods for modeling of contact interaction lead to significant errors in the obtained results or to overly cumbersome numerical models.Another important factor is the lack of adequate modeling of the contact conditions at the boundaries of the bodies.For instance, the impenetration condition for smooth bodies in linearized form is most commonly
At the same time, the surface rough layer has non-linear properties, which applies first of all to the "pressure-displacement" relationship.In addition, the contact may be as well established in the presence of a lubricant material or through a hydrodynamic layer in the interaction area.The machine parts surfaces can be subjected to chemical, mechanical or heat treatment as well.Elastic gaskets, films or seals can be placed between the bodies.This fact significantly affects the physical and mechanical characteristics of contact surface layers, which are based on micromechanical models, and can become considerably non-linear.The good example is the roughness of the machine parts surfaces.So far, the data of experimental studies or various simplified contact interaction models of the microasperities in the form of rods or semispheres are employed in order to derive these relations.In this case, the different models of local compliance, which are generated by the surface microasperities (or other sources), are established based on the obtained data.It should be noted that the models that describe the connection "normal displacement -contact pressure" are substantially non-linear.These models cannot be linearized without loss of physical adequacy and numerical accuracy.
Therefore, nonlinear components appear in the conditions of the contact interaction (impenetration).They represent physical and mechanical nonlinearity of the material characteristics of various layers, in particular, the roughness.As a result, the structural nonlinearity is supplemented by a physical one.The physical nonlinearity is present in the relations that represent the essence of the structural nonlinearity.It is necessary to develop new methods and approaches to solve such problems, which contain non-linear terms in the contact conditions.These circumstances form the actual scientific problem of developing the new algorithms for analysis of contact interaction of construction elements.They have to take into account not only structural, but also additional physical nonlinearity.

Literature review and problem statement
The analysis of the current state of contact research between complex-shaped bodies.As for today, contact interaction mechanics (or contact mechanics) has formed into a well-established research area [1].At the same time, numerous methods related to the discretization of the bodies by the finite element method (FEM) [2] and boundary element method (BEM) [3] have received essential development.The theoretical basis of these numerical methods for contact problems can be the variational formulations and boundary integral equations respectively [2].
As for variational foundations, they have received significant development on the basis of the variational inequality theory [2].Contact pressure and contact area are defined from the extremum conditions of some functionals.In particular, the Kalker's variational principle is one of the options of this type of formulations [5].This variational principle operates with a functional defined on the set of nonnegative distributions of contact pressure.In fact, the minimum of the functional is sought on this set.In addition, a variety of directions for the development of variational statements are possible [6].
The boundary element method has its advantages and disadvantages as one of the variants of discretization of boundary integral equations [1,3,4].On the one hand, the physical dimension of the numerical problem is reduced by one.On the other hand, it requires operations with densely populated instead of sparse matrices.
The question arises: how the physically nonlinear properties of the interface surface layer of the contacting bodies (in particular, related to their roughness) can be taken into account within the given mathematical model.It should be noted that several very diverse theories of roughness are available at this moment.
Contact micromechanics: approaches and models.The surfaces of real bodies at the microlevel are not perfectly smooth, but there is a random profile consisting of convexities and concavities.This means that the initial contact will not occur over the entire nominal surface, but only on its small part that will grow with the increasing force due to the microasperities deformations.
The theory of roughness and its influence on the nature of contact interaction is far from completion.The first attempts to create a model of rough contact are based on the application of Hertz's theory to individual asperity vertices.However, the obtained theoretical estimates did not coincide with the experimentally established law of proportionality between the force magnitude P and the contact area A. The contradiction was that not only the area of existing contact spots increases with the growth of the load, but new ones arise [7].To solve the problem, an approximation is proposed.According to it, the microasperities are located on larger asperities, and even smaller microasperities are located on these microasperities.Such model gave the desired proportionality in the limit value with increasing amount of the levels.The following development of the rough contact theory was aimed at obtaining not only qualitatively correct, but also quantitatively accurate estimates of behavior.The methods of statistical averaging were used for this purpose.A comprehensive theory that describes the statistical properties of the random surfaces is proposed in [8].However, even before its appearance, Greenwood and Williamson developed a contact model based on simplifying assumptions about the distribution of the asperities characteristics [9].In particular, the shape of the bodies was adopted as spherical, so that contact is established on circular areas, and the curvature is a constant for all peaks, regardless of height.Random distribution of the asperities curvature is taken into account to develop the model that is proposed in [10].Subsequently, individual statements of the model were revised in [11].It is proposed to use a two-parameter Weibull's distribution for describing the heights of the asperities [12].The contact of individual asperities was taken into account independently in the original Greenwood-Williamson theory.The combined influence of the deformations caused by the application of the contact forces at adjacent vertices was taken into account [13].Individual models contain plastic deformations [14][15].The presence of these deformations is usually excluded by the fact that the asperities do not increase the value of the maximum contact pressure relative to the mean value so much that it leads to significant plastic deformations.In addition, even if the plastic deformations appear, there is a reason to believe that all subsequent loads will be elastic after several loading cycles.The disadvantages of the models describing the behavior of interface layers given in [7][8][9][10][11][12] consist in the fact that they do not take into account the mutual influence of the adjacent microasperities deformation.Some progress in the development of the mutual influence of the microasperities is proposed in [13].
The further development of the rough contact theory is related to the fractal nature of the bodies geometry.The theoretical results [16] showed the significance of the scale-effect of the surface topography on the predicted contact behavior.Persson developed an in-depth model, in which substantiated the inverse-potential dependence of the contact pressure value on the value of the middle gap between the pressed surfaces.This behavior is expected when the surfaces are compressed moderately, and the approach leads to the development of contact on a large number of vertices, but far from the establishment of complete contact.The numerical simulations were carried out to verify the basic statements of this theory.This confirmed basic assessments qualitatively and generally also quantitatively.
The contact stiffness is another value that reflects the properties of rigid contact.In [17], a direct analogy between the problem of elastic contact and electrical conductivity is examined, it is directly related to the resistance of the contacting bodies [9,18].In accordance with the Greenwood-Williamson theory, as well as Persson's theory, the contact stiffness is directly proportional to the pressing force.However, a number of other studies indicate a power law with the exponent from 0.5 to one [19].The key parameter that affects its accurate value is the fractal dimension of the surface, which was expressed in an analytical evaluation [20].
The fractal nature of the geometry of real machine parts is taken into account in [16].However, the relations of the studies [13,16] apply only for moderate pressure.In order to develop the theoretical relations in [9,[17][18][19][20], the power-law dependencies between the displacements of the surface layers points and the contact pressure are proposed.A special case is the contact with roughness and adhesion, which is distinguished by a fundamentally different law of surface interaction [21,22].
The application of any aforementioned theories in the practice of engineering analysis requires significant changes of the simulation tools that are currently being used.FEM [1,2] and BEM [1,3,4] have limited possibilities for modelling non-linear properties of the interface layers.Kalker's variational principle [5] also needs to be developed and adapted for the case of nonlinearly elastic bodies.From the view of the variety of the non-linear models of roughness, it is desirable that a new formulation would be general, i.e. assume any arbitrary law of the dependence of the state variables within the interface layer.In [6], an additional term was proposed for the Kalker's variational functional, this term takes into account the linear compliance contribution of rough surfaces in general terms.Subsequently, this approach was applied to create the macromodels of the contact interaction and implemented using the author's methods for solving the system of contact equations and inequalities described in the papers [21,[23][24][25].
The analysis provides the basis for the following conclusion: a combination of the following several methods is best suited for the analysis of the stress-strain state of complex-shaped bodies with physically non-linear interface layers.The variational formulation of the Kalker's type principle, the boundary integral equation method and the discrete approximation of the unknown functions are used to investigate the contact interaction.These methods are presented as a partial sum of a series of basic functions with a local support.
Given the circumstances, it is considered expedient to develop existing author's methods and adapt them to solve the applied contact tasks [6].At the same time, their direct application is impossible as far as the system of solving relations is not only structurally, but also physically nonlinear.

The aim and objectives of the study
The aim of this research is to develop methods for contact analysis of complex-shaped bodies with account for non-linear interface layers.
The following objectives were set to achieve this aim: -to develop more advanced macromodels for contact interaction of complex-shaped bodies based on micromechanical approaches; -to develop additional gap method and additional compliance method to solve the obtained system of non-linear relations; -to study the features of the contact pressure distribution with the variation of the gap shape and the properties of the interface layer between the contacting bodies.

Model of contact interaction with account for interface layers with non-linear properties
Consider the model of contact interaction of complex-shaped bodies following [1,5].When studying the contact of smooth incompliant bodies in the absence of friction [1], the displacements of the surface points and the gap between them are considered only in the normal direction as the first approximation.In this case, the kinematic relations of the contact are determined.
In the general case, the accurate form of the gap between the bodies has to be taken into account.The coordinate system with its center (point O) traditionally located on the line of the action of the pressing force P (Fig. 1, a) is introduced for this purpose.The equation of each of the surfaces is ( , ), and the gap has the form of Both bodies are deformed and come into contact on a certain area in the deformed state under the action of the force P (Fig. 1, b).The approaches of bodies i d are not related to their deformation.The displacements i z u correspond to deformations caused by the action of the required contact pressure.The generally accepted way of setting the nonlinear relations for normal contact is as follows: The integral relation is known for an elastic halfspace.It establishes the relationship between pressure and normal displacements of the points of its boundary (Fig. 1, c) [1]: It can further be derived as: where , , are the Poisson's ratio and the modulus of elasticity for the material of each of the contacting bodies, while the pressure distribution ( ) , p ξ η and the contact area S are the sought-for unknowns.
The piecewise linear representation of the normal stress distribution is used.The unknown contact pressure function is approximated by the superposition of an array of pyramidal elementary distributions.Their vertices are located in the nodes of the regular mesh with the step с.This mesh consists of equilateral triangles.The unknown function is fully determined by a discrete set of nodal values of pressure n p : Fig. 2. Regarding the approximation of the contact pressure distribution: a -the regular triangular mesh and the pyramidal element of the pressure; b -the displacement distribution of points of the half-space surface under the pressure The following two approaches, the direct method (or the method of collocation) and the variational method (Kalker's principle) can be used to find the pressure values in the mesh nodes that best satisfy the boundary conditions [4].The application of each of the aforementioned methods requires the ability to calculate normal displacements (3) for the chosen pressure approximation.Obviously, this task is equivalent to the determination of the displacements z u corresponding to each of the basis loads.Due to the homogeneity, it is sufficient to calculate the values of surface displacements for the one unit pyramid with unit sides (Fig. 2) as where , is the hexagonal area of the pyramid element with a vertex at the node ( ) S is the hexagonal area with unit sides, (1)  p is the unit pyramidal distribution on it, and the "template" of the form of the displacement distribution for a unit pyramidal element (Fig. 2): ( ) , / d d .
A simple algorithm for calculating the values of this function in separate points as well as on the entire area [5] was proposed.The method that is based on the analytical determination of the subintegral expressions in (5) for each of the triangles was taken as a base.These triangles form the basis hexagon (painted in gray in Fig. 2).Thus, the accurate analytical relations are obtained for calculating the components of the "template" (5) (Fig. 3).This is a positive factor in terms of ensuring the accuracy of numerical modeling of contact interaction.Integer coordinate-subscripts ( , ) i j correspond to any node of the mesh, such that its radius vector is calculated as ( ) ~.
e e Also a sequential numbering of nodes can be introduced.The equality (4) can be rewritten for the nodal points using this indexing: ( ) The equation ( 6) reflects the calculation of the influence matrix of coefficients C, which relates the nodal values of displacements to the nodal values of contact pressure.
The direct solution method for non-Herzian normal contact problems for elastic smooth bodies consists in writing the contact conditions for a finite number of nodes.The obtained system of equations enables to find the nodal values of the contact pressure that satisfy them [5].When the same regular mesh is used for the pyramidal pressure elements and the collocation points, the contact conditions (1) can be expressed in the following discretized form: 0, node in contact; 0, node out of contact.The integral force identity is written in the form of: The system of equations and inequalities ( 7)-( 9) is solved by iterations.
Variational methods for studying the non-Herzian normal contact.The variational approach is based on the weak formulation of the problem.In the considered case, the Kalker's variation principle [4] is the most suitable framework.According to it, the full complementary work Ф reaches the minimum among all possible nonnegative distributions p for the real area and the contact pressure acting in contact bodies: where S is some large enough area of the half-space.The same basis functions as above (Fig. 2, a) are used for approximation of p.Using the quadrature formula ( ) it can be reduced to the next quadratic programming problem { } ( ) ( ) The expression (11) is the approximation of (10).
It is noteworthy that the nodal values { } 1 , Т n n p = obtained as a solution of (11), are identical to the result of applying the direct method and satisfy the conditions ( 7)- (9).The advantages of the variational approach are the availability of the conditions that determine the contact area and the contact pressure distribution simultaneously.Also, quadratic programming methods can be applied.The latter circumstance allows applying the variational formulation for the formal justification of using the collocation method.
The contact model of rough bodies.The model of an elastic layer was applied in the first approximation to take into account the influence of roughness on the contact interaction of complex-shaped bodies.The stiffness of the intermediate layer is equivalent to the properties of the rough surface layer in some sense.Then the considered body Ω consists of two parts: Ώ -the smooth elastic body and ∧ Ω -a rough layer that covers it (Fig. 3).
There is a dependence between the displacements of the points z u Σ of the surface Ś involved in the description of the contact interaction conditions: ґ .
Here, ґ z u is the displacement of a smooth body and is expressed through the integral relation (2), z u ∧ is the displacement of the rough layer described by the corresponding model.Various analytical models of the contact layer are considered in [18,19].Winkler's foundation [5] can be mentioned as the simplest possible variant: where λ is the compliance of the layer (or layers), which depends on the material properties and mechanical treatment of the surface layer.The model ( 12) is a rather simplified version of more adequate and accurate laws (for example, power-law dependency) between the deformation and the pressure in the rough layer.However, even this simplified form makes a qualitative change in the system of interacting bodies by introducing additional elements into it.The governing relations take the form of: 0, node in contact; 0, node out of contact, where nm nm nm C C Σ = + λd are the coefficients of the influence matrix, nm C are the coefficients of the compliance matrix, which are determined by the relations (6); λ is the total compliance of the rough layers; Adding diagonal terms to the complete compliance matrix does not affect its positive definiteness.The equation ( 9) is not altered.As a consequence, the resulting system of equations by its structure remains identical to the smooth body contact but only incorporates the contribution from the specially introduced Winkler's layer.
The presented relations are the starting point for the development of advanced models and methods for analyzing the contact interaction of complex-shaped bodies taking into account the nonlinear characteristics of the interface layer material.
The model with a nonlinear Winkler's layer between contacting bodies.The analysis of contact interaction of complex-shaped bodies with an intermediate nonlinear Winkler's layer requires further development of the previously proposed methods and models.As outlined those only considered the contact of the bodies that are either smooth or have roughness represented by just a linear elastic layer between them.The models proposed in the work above have a certain advantage, as they allow for gradual "build up" by introducing new additional factors right into the initial formulation.This is due to the fact that the relations adopted as a basis are nothing more than the compatibility conditions for the displacements of the points interacting when bodies get in contact.These geometric relations do not depend on the physical and mechanical properties of contacting bodies in the initial form and are strictly fulfilled in the actual state for any mechanical system.The influence of the elastic properties of the interacting bodies is represented in each case by the particular dependence of the surface points displacements on the contact pressure.This is the essential component that can be updated in the existing models.Meanwhile, the original principle and the model structure of the fundamental formulation are not changed in this case and remain essentially the same as for the smooth contact.The new components are simply added to the model, which introduces new physical effects but preserves the inherited initial mathematical structure of the contact problem.
Thus, not only a linear model, but also a micromechanical model of the general form ( ) ( ) can be used to represent the surface roughness or another elastic layer between the contacting bodies.
In particular, the power dependence of displacements on the contact pressure (for taking into account the properties of roughness) is described in a number of papers [1,5].Other relations can also be considered and implemented as long as they adequately describe the properties of the layers between the contacting bodies.
Then the equation system in the expanded and matrix form is

C p w p h n N c p P p Cp D p h p P p c
where ( ) D p is the vector of layer deformation components ( ); w p С is the influence matrix;
Thus, the initial system (14), which previously contained only linear components in the displacements compatibility conditions in the left part is updated with the nonlinear components ( ).
w p In other words, the structural nonlinearity is supplemented by a physical one.Because of this, it is not possible to carry out a direct linearization in the relations of the upper row (14) in the general case.This is a fundamental difference between this model (addressed as structurally and physically nonlinear in the work) from the traditional structurally nonlinear, but physically linear models.
The relations (14) or similar can also be obtained from the generalization of the Kalker's variational principle where the third component in the expression ( ) p Φ describes the deformation energy of a nonlinear elastic layer (it corresponds perfectly to the term  14) are obtained by using the same quadrature formula as in (11).Thus, it is possible to propose a universal method by which qualitatively new mathematical model can be easily derived for the contact interaction of the complex-shaped bodies system in the presence of nonlinear deformable interface.The proposed method consists in the formation of the functional ( ), p Φ that contains the energy of the elastic deformation of all system components.This procedure does not represent a significant complexity due to the fact that the functional ( ), p Φ is additive.At the same time, each additional component in ( ), p Φ can bring new types of nonlinearity.This in turn naturally resolves all the issues related to the existence and uniqueness of the solutions to the problem, as well the convergence and accuracy of the developed discretized models.

Solution methods for structural and physically non-linear problems of contact interaction
The mathematical model of contact interaction of smooth or rough complex-shaped bodies and of the bodies with nonlinear elastic layers between them was described above.The structurally-physically nonlinear system of relations was obtained as a result.It is necessary to develop new methods for solving this system of relations, since the traditional ones are only applicable to the solution of contact problems with linearly elastic components.
As an example, a nonlinear model of an elastic layer can be accepted as , where , s λ are some empirical or computable parameters.The compatibility relations can be formally reduced to a system of nonlinear equations as a result.Thus, the problem of developing methods for solving such tasks arises.
The augmented gap method.If the displacements compatibility relations that enter the system ( 14) are written as a subsystem , ,..., , then formally it can be represented as

( ). h h p h q h h p
Then, the relations (18) formally repeat the linear relations in (20), but its initial gap h is augmented by some compensating components ( ). h p D These components can be called the additional or augmented gaps by analogy with the additional load method.
The equation ( 18) is a nonlinear operator equation for which the real distribution of the contact pressure is a fixed point of the operator of the complete equations system under conditions 0. p ≥ Thus, the iterative process can be organized as The variables Р, q, h, d are refined iteratively at the final step (^^^) in (19).Herewith, (0)  h is the vector of initial (nominal) gaps between the surfaces of the contacting bodies.Iterations are carried out within the cycle (***) (:) (***).
→ → The condition for the completion of this process can be either a criterion in terms of the increments of additional gaps, or the pressure: / , where ( )  Thus, the solution of the problem is equivalent to the solution of the problem (7) for the contact of the smooth bodies, but with the corrected distribution of the gaps.
The method of variable coefficients of influence.Presenting the displacement compatibility equation in the system (14) in the presence of the Winkler's layer on the boundary of bodies ( 16) in the form of ( ) and introducing the notation = λ one can rewrite (by analogy with the method of variable parameters of elasticity) the equation ( 21) in the form of An iterative process can be arranged to solve (22): In (23), the iterative process goes in the sequence (***) (:) (***).

→ →
The termination of the iterative process can be implemented on the condition of stabilization of the iterative approximation of the set of nodes in which the contact is realized: ( ) ( 1) .J J β β− = The presented algorithm (23) differs from the algorithm ( 19) by the fact that not the right-hand side of the system of equations with invariable components of the submatrix С is being changed, but on the contrary the matrix itself is changed while the right-hand side is constant.
The proposed approach can be interpreted as follows: such a distribution on the nodes of the compliance coefficients mesh ( ) n n J λ ∈ of the linearly elastic Winkler's layer with an irregular stiffness on the contact area is determined in the solution (22) by the algorithm (23) that the solution of the problem with the matrix : , gives the same distribution of the contact pressure, as in the case of a nonlinearly elastic layer (16).
In general, the developed methods for solving problems of the contact interaction analysis of the complex-shaped bodies with a nonlinearly elastic interface layer are adapted to solve the obtained systems of equation and inequalities.This determines the advantage of the developed methods compared to the known traditional methods for solving nonlinear equations.

Test problems solution and results analysis
The methods and models of contact interaction analysis developed for complex-shaped bodies in the presence of a nonlinear elastic interface layer need to be tested by solving a set of representative problems.In particular, the convergence of the numerical solution using the developed iterative procedures has to be evaluated.The test problems chosen for this purpose involve contact of two bodies limited by paraboloids with nonlinear elastic layers between them (Fig. 4, a).Varied are the elastic characteristics of the nonlinear layer in terms of the function ( ).
w p In particular, the following types of nonlinearities were considered (Fig. 4, b-d where max p λ is the maximum contact pressure in the distribution obtained during the solution of contact problems without limitations cr p (the case 0 a α = corresponds to the contact of smooth bodies).
As it can be seen from the presented distributions, there is a gradual transformation of their shape with the variation of the coefficient a α in the interval [0; 1] (Fig. 6, k is the iteration number).The zone of transition between the regions І and ІІ is clearly visible (Fig. 7).There are different distributions of the contact pressure q in regions І and ІІ in the zone of 0,5, a α ≈ and the transition zone is between them.
The distributions q are close to the limit at approximation a α to 0 or 1 (respectively, for 0, The change of the relative maximum contact pressure during the iteration process (Fig. 4, b) is illustrated in Fig. 8.It is seen that the iterative process converges sufficiently quickly.In addition, the contact pressure level varies smoothly, but within a sufficiently wide range when varying .The contact pressure distributions in the variation of the coefficient p The contact pressure distributions obtained not only by the augmented gap method, but also by the method of variable compliance parameters are shown in Fig. 6.It is seen that the results obtained by different methods generally coincide.
Thus, as it can be seen from the obtained results, the coefficients , a α , b α c α strongly affect the character and the level of contact pressure distribution.As a result, the stressstrain state of contacting bodies can be significantly affected by varying these coefficients.
The solution of the non-linear equations systems (13) and inequalities systems (14), which compose the discrete contact conditions, was carried out by using a software developed in previous studies, in particular [19].
The numerical algorithm and its implementation have tested on a number of test tasks, in by comparing them with known analytical theories and the results of finite-element modeling in ANSYS.
In addition, the error behaviour of the approximate solution is studied.It displays constant and rapid convergence, and hence the accuracy of the presented results.

Discussion of the proposed model capabilities
The problem of contact interaction of linearly elastic bodies with a nonlinearly elastic layer between them has been reduced to a structurally and physically nonlinear problem by means of the boundary integral equation method and Kalker's variational principle.The obtained mathematical model is related to the system of equations and inequations in contrast to the conventional formations.These equations and inequations make up the conditions for compatibility of displacements, which contain nonlinear components with respect to the non-negative contact pressure.Thus, the structural nonlinearity, mathematically represented in relations of the inequations type, is supplemented by the physical one, which gives rise to the appearance of nonlinear components in the equations of compatibility.It complicates the mathematical model, but makes it more adequate and describes the behavior of real objects more accurately.On the other hand, it requires the development of new methods for solving the obtained system of relations.
An alternative way of creating a mathematical model of a contact interaction of systems of bodies and layers is to use a variational formulation.In this case, in addition to components of the initial level (the energy of the elastic bodies deformation), the components that describe the contribution of other components can gradually be added.Afterwards, a direct procedure for minimizing the functional is applied, or a system of equations and inequations is obtained from the minimum conditions and constraints respectively.
For the solution of the formed structural-physically nonlinear relations, the augmented gap method has been developed.It reduces the initial structurally and physically nonlinear problem to a sequence of structurally nonlinear but physically linear problems, each having a particularly evaluated distribution of the fictitious supplementary gap.Since effective iterative algorithms for solving such tasks were developed previously, this allows solving effectively and accurately physically nonlinear problems.
In addition to the augmented gap method, the method of variable compliance parameters has been developed.It consists in reducing the initial physically nonlinear problem to a sequence of physically linear ones through the specially introduced stiffness of an equivalent linear Winkler's foundation that is distributed not uniformly over the contact spot.
Compared to the case of smooth bodies, the additional compliance of rough surfaces leads to an increase in the contact area and a decrease in the maximum value of the contact pressure.At the same time, the considered non-linear laws of the pressure dependences (in particular bilinear and power ones) and an approach in the layer are more physically adequate in view of the finite ability of asperities to the compression.

Conclusions
1.The new approach for creating mathematical models of the contact interaction of elastic bodies with nonlinearly elastic layers between them, as well as taking into account other factors, is proposed.It is based on the gradual enhancement of the structure in accordance with the introduction of any new particular factor.The initial core in this sequence is the displacements compatibility equation for surfaces points of smooth bodies (in the actual state) in terms of contact pressure.Further, the components representing an additional factor are added to this core.As a result, the system of relations that describes the influence of various factors in total is obtained.
2. The approach based on using the boundary element approximations of the sought-for contact pressure on a regular triangular mesh is used to discretize the mathematical model.As a result, the problem is reduced to minimizing the convex function on a convex set of nodal values of contact pressure.
3. New regularities in the contact pressure distribution between the complex-shaped bodies for different dependencies ( ) w p have been obtained.It is determined that the obtained pressure distributions have two characteristic regions with a transition zone for the case of an "elastic-rigid" characteristic of the properties of an interface elastic layer.At the same time, changing the "threshold" of the transition from the elastic to the rigid branch of this characteristic leads to a continuous change of the obtained contact pressure distribution.It is established that the obtained pressure distributions have a more smooth transition to zero distribution outside the actual contact area compared with the case of a linear elastic layer for the problem with the "root" characteristic that describes the properties of an interface elastic layer.

Fig. 1 .
Fig. 1.Regarding the normal contact of complex-shaped bodies: a -the representation of the local gap between the contact bodies; b -the deformation of the bodies and the formation of the contact area under the action of normal forces; c -the displacement of the half-space under the action of normal forces u x y h x y S x y nd S x y of the initial gap.Another condition is the nonnegativity of the pressure that restricts the unknown nodal values

Fig. 3 .
Fig. 3.The model of the rough complex-shaped body: S ¢ is the boundary of the smooth body ;Ω¢ S ∧ is the border of the complex-shaped body case of a linearly elastic layer).The relations anal-ogous to (

D
respectively, after which the mentioned iterative process stops; • is some vector norm.

Fig. 4 . 3 Ω
Fig. 4. Regarding the test problem: a -the scheme of interaction of paraboloids 1 2 , Ω Ω with a nonlinear layer 3 Ω placed between them; b -the elastic-rigid dependence ( ); w p с -the "root" dependence of the form ; w p = λ d -the bilinear dependence ( ) w p The distributions of contact pressure depending on the threshold pressure value at which the transition to the horizontal section of the dependence ( ) w p (Fig. 4, b) takes place are presented in Fig. 5.The threshold is set by the value of the coefficient cr max /

aα
The changes of the contact pressure distribution with varying parameter a α are shown in Fig. 9.The evolution of the contact pressure distribution is noticeable, as well as the transition zones between branches І and ІІ the dependence ( ) w p (Fig. 4, b).The contact pressure distributions when root dependence ( ) w p (Fig. 4, c) are shown in Fig. 10, 11.The coefficient max max / b p p λ α =is being supervised, where max p is the maximum contact pressure at the "root" distribution ( ), w max p λ is the same one, but for the linear Winkler layer.It is seen that the "root" model (Fig.4, b) generates a more smooth contact pressure distribution, especially in the peripheral zone, in which the pressure level decreases sharply.

Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .
Fig. 8.The change of the relative value of the maximum contact pressure during the iteration process Fig.4, d) are also shown in Fig.5.As follows from the analysis of the obtained distributions, the contact pressure distribution tends to be the same as for the elastic-rigid model (Fig.4, b) when c α approaches to zero.The contact pressure distributions have a lower level in the middle part and wider contact areas accordingly when c α approaches to / 2.