A NUMERICAL METHOD FOR AXISYMMETRIC ADHESIVE CONTACT BASED ON KALKER’S VARIATIONAL PRINCIPLE

The phenomenon of adhesion occurs due to the action of molecular forces when the surfaces of solids come close within their range. As a result, not only compressive but also tensile tractions can be transmitted through the contact interface. The surface potential that favors expansion of the contact area can drive additional deformations required to close the gap between the bodies. As a consequence, even when the loading is removed a remnant contact is established and some additional pulling force is required in order to disengage the bodies. The classical Johnson-Kendall-Roberts theory [1] for the adhesive contact of parabolic bodies determines the actual contact radius by the minimum of full internal energy that includes the energy of elastic deformations and the energy of free surface loss. Similar analytical solutions are possible for cases with different geometries such as the contact of a paraboloid with an axisymmetric wavy surface [2]. Another well-known approach suggested by Derjaguin, Muller and Toporov [3] assumes the occurrence of attraction in a vicinity just beyond the contact area where the separation between bodies remains within a certain threshold. The advantage of this model is that it accounts explicitly for the finite range of molecular forces. On the other hand, it entirely disregards the contribution of these tractions to the elastic deformations of the bodies, which is a definite drawback. In general, analytical models are restricted to certain geometry and are frequently based on approximate representation of the nature of adhesion. Hence, numerical methods of analysis are required to predict adhesive behavior in case of arbitrary geometry or more sophisticated physical models of surface interaction.


Introduction
The phenomenon of adhesion occurs due to the action of molecular forces when the surfaces of solids come close within their range.As a result, not only compressive but also tensile tractions can be transmitted through the contact interface.The surface potential that favors expansion of the contact area can drive additional deformations required to close the gap between the bodies.As a consequence, even when the loading is removed a remnant contact is established and some additional pulling force is required in order to disengage the bodies.
The classical Johnson-Kendall-Roberts theory [1] for the adhesive contact of parabolic bodies determines the actual contact radius by the minimum of full internal energy that includes the energy of elastic deformations and the energy of free surface loss.Similar analytical solutions are possible for cases with different geometries such as the contact of a paraboloid with an axisymmetric wavy surface [2].
Another well-known approach suggested by Derjaguin, Muller and Toporov [3] assumes the occurrence of attraction in a vicinity just beyond the contact area where the separation between bodies remains within a certain threshold.The advantage of this model is that it accounts explicitly for the finite range of molecular forces.On the other hand, it entirely disregards the contribution of these tractions to the elastic deformations of the bodies, which is a definite drawback.
In general, analytical models are restricted to certain geometry and are frequently based on approximate representation of the nature of adhesion.Hence, numerical methods of analysis are required to predict adhesive behavior in case of arbitrary geometry or more sophisticated physical models of surface interaction.

Literature review and problem statement
The most fundamental approach consists in representing the microscopic structure of the bodies and modeling both the elastic behavior and contact interaction by means of molecular mechanics [4].This requires determining atomic forces between all the particles in the system, which has very high computational cost.
A method of intermolecular force homogenization allows upscaling the problem to the continuum level and perform effectively coarse-grained finite element analysis of adhesive contact [5,6].Molecular potentials can be also incorporated into the boundary element method, where they constitute
A heuristic criterion of surface detachment in adhesive contact is proposed in [10,11].According to it, contact is broken abruptly at a discrete element of a rectangular mesh when a certain level of tensile stress is achieved.The critical value is determined by the specific energy of adhesion, elastic material properties and the mesh size.However, this method also has high requirements regarding the mesh size for which accurate results are obtained.Moreover, the proposed algorithm does not ensure the uniqueness of the solution.
When contacting surfaces are rough, the interaction is localized initially around the asperities.The entire distribution of the contact tractions can be approximated then by a discrete system of localized contacts which is the main idea of the so-called multiasperity methods [12][13][14].On the plus side, it enables analysis of essentially large systems with relatively few unknowns.However, the discrepancy caused by simplifying assumptions is hard to control.In particular, it increases once the distance between separate contact spots shrinks so that they actually merge together.
It can be concluded that an approximate method for analysis of surface adhesion with lesser computational cost and acceptable accuracy is still required.It is proposed to develop it based on the existing boundary element methods, using in particular generalized minimum principle for adhesive contact and its formulations proposed in [15].This requires adaptation of the variational statement for the numerical application in the general case of axisymmetric contact.As long as the developed method is meant for modification of existing boundary element programs, Abel transformation is avoided unlike [16].It will be used for analysis of adhesive properties of various bodies with surface geometry met in practical applications.In particular, the influence of roughness on the strength of adhesion is a crucial question.

The aim and objectives of the study
The aim of this work is an effective and versatile realization of the boundary element method for adhesive contact problems.It is meant to preserve thermodynamical consistency of fundamental Johnson-Kendall-Roberts model.It suggests generalization of the variational formulation by J. Kalker which adds a separate energy term representing the contribution of adhesive interaction.This in turn allows for a standard approach to approximate solution by the Ritz method.
To achieve the set aim, the following tasks need to be solved: -formulate a variational principle of stationarity of total complementary energy in surface form with respect to the unknown contact domain and the distribution of contact pressure within it; -build a discretization of the variational problem for the case of axisymmetric contact with simply connected circular contact spot and implement a numerical method to solve the system of nonlinear equations; -perform analysis of model problems, estimate the accuracy of the method and verify the obtained results.

Variational formulation and numerical solution of the adhesive contact problem
The following mini-max principle is proposed for the problem of contact with adhesion max min [ , ], where C -the sought-for contact domain, p -the unknown contact pressure, u[p] -the total displacements at the boundary of the elastic semispace that are evaluated in the form of integral with the singular kernel derived from the fundamental solution of the Boussinesq problem where E * -the equivalent elastic modulus for both bodies, h -the initial gap, δ -the given relative approach of the bodies, γ -the specific surface energy.
The internal minimization with respect to the contact pressure p with the contact domain being held fixed involves only the first two terms of the functional (2) which compose the expression of the total complementary energy in the original Kalker principle [17].However, the values of the contact pressure in (1) contrary to the above-mentioned principle that is valid for unilateral contact are not limited to strictly positive values.Naturally, adhesion implies that the bodies may attract to each other over some part of the contact surface.Thus, the variation of the function p is arbitrary and extremum conditions come as They are satisfied as long as the following equation for the elastic displacement, the initial gap and the relative approach hold in the entire domain of the variable function which totally corresponds to the definition of integration domain C as the contact area.It's easy to deduce that the result of inner minimization of the complementary energy e Φ is the opposite to the energy of elastic deformations: The contribution of the adhesive forces enters the functional (2) similarly with the opposite sign.When contact is established over some part of the surface C, it eliminates the free surface and the corresponding energy loss is d .
C S -g ∫ Accordingly, the problem stated in ( 1) is in total agreement with the variational principle formulated by Johnson [18] and Kalker [15], namely: the true contact area and the true contact pressure distribution that satisfy contact conditions deliver minimum to the total energy of the system.The major obstacle to the practical application of this principle consists in mathematical formalization of the variation of the contact set C in which the unknown and variable contact pressure function p is defined.Nonetheless, there are several possible ways to build approximate methods for adhesive contact based on the problem statement (1).
This requires first to perform approximation of the involved functions, namely the contact pressure and the initial gap distribution.Each can make use of an individual basis in the functional spaces: It should be noted that the exact functions together with the respective approximations are bound to the domain C. In case when the utilized bases functions are built on a mesh, the discretization will need to assign appropriate locations to the corresponding nodes x a and x j in accordance with the shape of the contact domain.The ultimate discretized version of the functional (2) will have the following form: where p a and h j are the nodal values of the variable contact pressure and the known initial gap between the bodies, whilst the bilinear form coefficients are computed for the given choice of the basis functions as The values of the coefficients B ab , W aj of the functional (9) together with the gap h j depend on the contact area.When it changes, the nodal locations x a and x j are changed as well, which is not trivial in case of arbitrary shape of contact domain.One can avoid these difficulties in case of axisymmetric contact for which a simple and effective numerical method can be proposed.When the contact domain is restricted to the circular area aligned with the axis of radial symmetry, its change is determined by a single numerical parameter s: Due to the symmetry, all featured functions can be defined on the one-dimensional radial axis r x = : ( ) ( ), p x p r = ( ) ( ). h x h r = Introducing the affine transformation r sr = to the dimensionless coordinates r assigned to the interval [0,1], one can resolve the problem with the locations of nodes of discretization, the appropriate choice of the basis functions and the evaluation of the approximate functional.The proposed mesh is defined by a linear scaling on the segment [0,1].Piecewise constant basis functions for the contact pressure and piecewise linear approximation for the gap are identically defined in actual and dimensionless coordinates Such approach allows obtaining the final expression for the discretized functional and computing its coefficients as 3 , The consuming computations of dimensionless , ab B aj W need only to be accomplished once.To a great convenience and economy, the actual values of B ab , W aj are updated by simple multiplication at any subsequent solution step.This advantage leads to substantial acceleration of the numerical procedure.
Stationarity conditions of the approximate functional (16) form a system of nonlinear equations which define the nodal values of the contact pressures p a and radius s of the circular contact spot.This system is solved numerically by the Newton-Raphson method.

1. Contact of an elastic sphere with a plane
Contact of a sphere with a plane or equivalently of two spheres is well described by the Hertz theory [18].The axisymmetric initial gap in this case is approximated by a quadratic function 2 ( ) , 2 where R is the combined radius of curvature.The analytical expression of the ellipsoidal distribution of the contact pressure is obtained from the exact solution of the appropriate singular integral equation.The actual size of the contact spot is determined from impenetration condition that needs to be fulfilled outside.According to this theory, the pressing force, the relative approach and contact radius are related as: , 16 Unlike the Hertz theory, the proposed variational formulation does not account explicitly for the displacements beyond the defined contact area.However, the stationarity conditions for the functional (2) as proved in the work [15] provide the expected behavior of the approximate solution.
In particular, the contact pressure values remain positive and converge down to zero on the edge of the contact zone as can be seen in Fig. 1.Comparison of the obtained numerical results with the theoretical estimate of the compressive force and the contact area size ( 21) is presented in Fig. 2.Even a discretization with a few nodes is sufficient for gaining acceptable tolerance.Nevertheless, certain discrepancy in the computed contact area can be noticed.It is explained by the low sensitivity of the approximate functional ( 16) and its derivative with respect to the variable s when there's no adhesion.The Johnson-Kendall-Roberts (JKR) model for axisymmetric contact of a sphere in the Hertzian parabolic approximation has the variational basis [1] that is essentially similar to the principle used in this work.It differs though from the proposed approximate solution method as long as in case of the quadratic initial gap (20) the contact pressure distribution is known and expressed analytically for the fixed approach and the size of the circular contact domain.
The following exact relations are therefore established by the JKR model It can be easily noted that this solution reproduces the Hertz theory formulas (18) once 0, g = which means that there's no surface adhesion.Otherwise, the contact behaves completely different.First of all, the contact pressure distribution is not exclusively positive in the entire interaction domain anymore.Moreover, the negative values grow indefinitely with approach to the edge of the contact spot as can be seen in Fig. 3. Secondly, the relations (19) between the main quantities unlike in the Hertzian case lack strict monotony.Stable contact spot between the two bodies pressed together with a sufficiently large force narrows down when they get pulled apart.The force of contact interaction turns from positive to negative and attains its minimal value of  4 .
In a force controlled experiment, this is a moment when the bodies will undergo abrupt and complete detachment.In case of kinematic control, the distance between the bodies can be further increased and stable contact will be maintained up to the moment when the contact radius drops to 2/3 / 3 c s [18].The attractive force will become 5 9.
c P -at this point.Computation of the entire loading curve including the unstable portion requires numerical continuation.A spherical arc-length control is proposed to solve numerically the nonlinear stationarity conditions (19).It includes increments of the approach, the total contact force and the size of the contact.The latter parameter is essential for the robust convergence of the iterative procedure.

2. Contact of an elastic sphere with an axisymmetric wavy surface
Qualitatively different adhesive properties can be observed for more complex contact geometry.Roughness is known to have a drastic effect on the strength of adhesion [19][20][21][22].A model problem that allows studying this phenomenon is proposed in [23].
It considers contact of an elastic sphere with a wavy surface of convolution as shown in Fig. 6.As long as the radial symmetry is preserved, the problem may be solved by the proposed numerical method as well.The input function of the initial gap for this case takes the form where A is the amplitude of the cosinusoidal profile of the wavy surface, λ is the corresponding wavelength.Under certain conditions, the derivative ( ) h r ′ will stay positive within the entire body interface.For such geometry, unilat-eral contact without any adhesion takes place over a simply connected circular domain.A sufficient condition can be established as the inequality ( ) which satisfy this inequality is shown in Fig. 7.As can be seen, the contact tractions remain positive everywhere, thus no interruption of contact will occur within the circular contact zone.It should be noted that both the width of the contact area and the distribution of the contact pressure differ from the Hertzian solution.The waviness causes stress concentrations around the asperities or rather the circular ridges of the cosinusoidal profile and accordingly traction depressions in between.Obviously, with the increase of the amplitude A, this deviation from the Hertzian case might get such magnitude that the contact pressure will cross the boundary of positive values over a certain fraction of the contact which will be automatically broken there.The reduction of the wavelength λ will have a similar effect.Nevertheless, increasing sufficiently the pressing force will always eliminate these gaps and establish full contact.For any present roughness, the contribution of the sphere curvature to the slope of the gap will become dominant once the contact is expanded at sufficient distance and will ensure the inequality ( ) 0 h r > ′ over the major part of the contact surface.This statement is supported by the results presented in Fig. 7, b that were obtained numerically for 0.1 A l = and 0.1 l = R in the assumption of contact domain connectivity.The above-presented analysis may cast an impression that surface roughness exclusively impedes contact.However, its effect on the adhesive properties as outlined in the work [2] is overwhelmingly favorable.Consider similar to this study the case of 0.01, for the adhesion parameter value ( ) The computed equilibrium curves for the contact force and the size of contact are shown in Fig. 8. Since the introduced roughness is relatively small, the deviation from the JKR model which covers the contact of perfectly smooth bodies is also not too significant.From a state of pressing characterized by conditions 0, d > 0, P > chosen as a starting point, the separation takes place gradually with decreasing approach δ and contact force P.This process may be driven by a controlled change of one of these two parameters in which the other case is the reaction.If the force is the controlled variable than the abrupt separation will occur at a point A where the value of P attains its absolute minimum.In the case of displacement driven detachment, the contact will be broken at a different point B on the equilibrium curve with the maximal distance -d.It can be seen that in both regimes a larger tensile force or accordingly a larger separation distance are required to detach bodies from each other compared to the case without roughness.This effect is even more notable when the roughness amplitude is amplified.The corresponding example for 0.05, A l = 0.05 R l = and 0.025 g = ′ is presented in Fig. 9.In this case, the increase in the force of separation is more than 50 %.Besides the load-separation curve becomes distinctively non-monotone.As a result, the detachment will take place in an unstable manner accompanied by a series of abrupt jumps of state of contact.This introduces a qualitatively new dissipative mechanism that is not present for smooth contact in the original JKR model.The described behavior is illustrated in Fig. 10 for both types of loading.In particular, for the displacement control the unstable jumps will occur at every point of the equilibrium curve where it is tangent to the vertical lines const d = (Fig. 10, a).During each such event, the contact spot will shrink instantly by a finite amount.The contact force will accordingly drop to the next stable value on the curve.The state of contact will alter in a similar manner along the horizontal lines const F = (Fig. 10, b) in case of the force controlled process.These critical events can be clearly pinpointed to the moments when the edge of the circular contact passes through the asperities of wavy profile.The next stable location of the boundary after the instant jump will be found on the following cosinusoid wave at a distance of order the wavelength λ.There is an irreversible loss of elastic energy which is freed upon every instability and abrupt change of state.Thus, the roughness is capable to increase the toughness of adhesive contact by this dissipative mechanism.

Discussion of results obtained with the proposed numerical method
The proposed variational formulation is based on the Kalker's variational principle of stationarity of total complementary energy.Compared to the general case, the variable contact domain can be determined without difficulties when the geometry is axisymmetric.This allowed constructing a simple and effective approximate solution method.The proposed numerical implementation benefits from the one-time computation of all quadrature coefficients that are updated with the variation of contact radius at each new iteration by simple scaling transformation.The arc-length method with control that includes increments of both contact approach as well as the contact radius is employed to compute the unstable sections of the loading curve.The advantage of the method is that it predicts continuous evolution of the contact domain together with its area.This allows employing it to study the effect of surface roughness on the contact interaction.In particular, the continuous solution provides derivatives of such quantities as the contact area and the relative approach of bodies which is essential for the estimate of electrical and thermal conductivity as well as contact stiffness.
In should be noted that for certain values of roughness and adhesion parameters, the contact may initially be established over a series of separate circular rings, which violates the basic assumptions.Nevertheless, further compression of the bodies to each other will inevitably force these rings to merge into a connected area of full contact.Hence, leaving the process of approach beyond the scope of analysis, it is still possible to use the proposed method to study their detachment.
The further development of the proposed approach is sought in extension to the case of disconnected contact and contact geometry different from axisymmetric.This requires variational apparatus that is capable of handling arbitrary two-dimensional contact domains.A possible solution can be obtained by level-set formalism previously applied to contact of a membrane with a rigid obstacle [24].The method can also be applied to other types of surface interaction different from the molecular adhesion, such as gas pressure and liquid bridges on the interface [25,26].

Conclusions
1.The work proposes a variational formulation of the axisymmetric contact problem for bodies with dimensions essentially larger than the contact zone with an account for adhesion.In this statement, the contact pressure and the radius of the circular contact domain are the unknown variables.They are determined by the principle of the minimum of total complementary energy derived from the well-known Kalker's principle for unilateral contact.The problems related to variation of contact spot are resolved due to the restriction to the axisymmetric case.The circular shape is determined by a single scalar parameter, while the domain of the sought-for contact pressure distribution function is defined by a simple affine transformation.
2. The boundary element method is used to obtain the approximate solution of the posed problem.The system of nonlinear equations determining the nodal values of the contact pressure and the contact radius is derived from the stationarity conditions of the discretized functional of total complementary energy.A numerical procedure to solve these equations by the Newton-Raphson method is developed and implemented.Exact derivatives of the left-hand side of equations and the spherical control are obtained and used for this purpose.The numerical procedure demonstrates stable and rapid convergence.Even when the initial approximation of the dimension of the contact spot has a threefold excess over the true value, the method converges just in 5-6 iterations.
3. The effectiveness and accuracy of the proposed method are verified by solving a series of classical problems.The method exhibits linear convergence to exact solutions of the Hertzian and Johnson-Kendall-Roberts problems with just a couple of dozen of nodes required to reach high accuracy.The complete response curves for the contact of an elastic body with a wavy surface are computed including the unstable sections are computed with the help of the implemented numerical routine of arc-length control.The obtained results complied for this case qualitatively and quantitatively with the analytical Guduru theory.Thus, the proposed numerical method was proved capable to compute contact behavior of axisymmetric bodies of complex shape in the presence of adhesion.

Introduction
The present-day production of construction machines is impossible without conducting a large test cycle before their delivery to the customer.Tests start from prototype models and end with serial machines.Prototype tests play an important role in improving the construction machine design as they can reveal drawbacks and improve reliability of the machine design.

SUBSTANTIATION OF ADEQUACY OF LOADING CONDITIONS AT BENCH AND FIELD TESTS OF CONSTRUCTION MACHINES L . P e l e v i n
PhD, Professor* Е-mail: pelevin_leonid@ukr.net

Y e . G o r b a t y u k
PhD, Associate Professor* Е-mail: gek_gor@i.ua

G . M a c h i s h i n
PhD, Associate Professor* Е-mail: ma4ichin@ukr.net*Department of Construction Machinery Kyiv National University of Construction and Architecture Povitroflotsky ave., 31, Kyiv, Ukraine, 03037

Fig. 3 .Fig. 5 .
Fig. 3. Computed contact pressure distribution compared to the Johnson-Kendall-Roberts modelThe adhesion affects the convergence of the approximate solution with the increase in the number of discretization nodes.The exact solution has a square root singularity at the boundary of the contact pressure domain.Compared to it, the Hertzian distribution is finite-valued, although its derivative has an analogous singularity of the same order.Both cases have different tolerance of approximation by piecewise constant functions.The resulting convergence of the proposed method is presented in Fig.5.It can be concluded that