1 Introduction
The analysis of the tangential tractions arising during the frictional contacts between a rigid and an elastic body plays a central role in physics and engineering. Mathematical models based on different friction laws have been developed during the past years in parallel with the progresses in the field of contact mechanics barber2020 ; vakis . In general, due to the intrinsic characteristics of the frictional contact problem, non uniqueness issues might arise, depending on the coefficient of friction and the regularisation of the Coulomb law adopted klarbring1990529 ; raous1 ; raous2 ; raous3 .
A rigorous analytical treatment of this problem permits to obtain a solution only in few selected cases, characterised by simple contacting geometries and constitutive behaviours. For example, starting from the general class of problems addressed in the first paragraph, if the attention is restricted to , small displacements and the deformable material is elastic and isotropic, the formal problem results in a set of two integral equations, coupled by the first and second Dundurs’ bimaterial constants and , equipped with a proper friction law barber_contact . If the attention is further restricted to instances in which the halfplane hypothesis can be invoked, the problem reduces to two integral equations coupled by the parameter only. Under these assumptions, the effect of coupling causes a vertical pressure to induce also horizontal displacements, thus influencing the tangential tractions distribution which in its turn modifies the vertical displacements, having so an appreciable effect on the vertical tractions and the contact area distribution. The halfplane approximation is particularly relevant since, under this assumption, problems involving also the contact of two elastic bodies can be addressed, via the introduction of adhoc elastic moduli.
Under these hypotheses, the pioneering investigations of Cattaneo cattaneo and Mindlin mindlin are particularly relevant. These studies seminally concern with the solution of the problem of two contacting elastic spheres. They independently showed that the tangential tractions distribution caused by friction could be expressed as a superposition of two normal tractions distributions. Later, still within the framework of uncoupled approaches, Jäger jager1 ; jager2 and Ciavarella ciava1 ; ciava2 extended CattaneoMindlin results to the contact problem of halfplanes with generic non compact boundaries under oblique loading. In goodman1962 , Goodman introduced the hypothesis of semicoupling, taking into account only the effect of normal tractions on horizontal displacements, but not viceversa, allowing the independent treatment of the vertical pressures with the usual methods employed for frictionless problems, and then solving separately for tangential tractions. This framework has been exploited by Nowell et al. nowell who obtained the analytical solution for the special case of a quadratic indenter acting upon an elastic halfspace. In the same reference, they also developed a numerical scheme based on the inversion of the governing integral equations for solving the related fully coupled problem with generic tangential loading histories. Under the hypothesis of full coupling, Spence spence solved the problem of the indentation of an elastic halfspace by an axisymmetric punch for a monotonic normal load. In particular, he solved the problem for a flat punch and a Hertzian indenter, and proved, thanks to a property of the governing equations, that the solution for a generic power law profile could be directly derived from that corresponding to the flat punch. As a consequence, under the same intensity of the normal load, the extension of the slip area is found to be the same for every indenter with a powerlaw profile.
More recently, numerical methods have experienced a considerable development in order to overcome the limitations of the analytical approaches. Without the aim of providing an exhausting literature review on this topic, credit to the main contributions in the development of this subject can be given to Klarbring klarbring1 and Klarbring and Björkman klarbring2 , Kalker kalker1 ; kalker2 and Ahn and Barber ahn . Numerical schemes based on the Boundary Element Method (BEM) have been proposed and widely used for the quantitative evaluation of tangential contact problems under arbitrarily complex contacting geometries. This framework has proved to be very efficient from a computational point of view, since it only requires the discretisation of the domain boundary, without accounting for the bulk. Many implementations are available, see e.g Zhao et al. zhao , Pohrt and Li pohrt and Willner willner . On the other hand, BEMbased algorithms are limited by the assumptions of linear elasticity, halfspace approximation and homogeneity of the materials. Their extension to inhomogeneities nelias , nonlinear interface constitutive response Rey2017 ; Popov2017 ; Popov or finite size geometries finite ; vollebregt2 are sometimes possible but yet exceptions.
The Finite Element Method (FEM) represents today a suitable modeling tool for overcoming the inherent shortcomings of BEM, since it can easily account, for example, for geometrical or material nonlinearities and finitesize geometries HPMR ; PEI . A shortcoming of FEM applied to both frictionless and frictional contact problems is represented by the treatment of complex contacting geometries, since a discretisation of both the interface and the bulk is requested by the method. In case of very irregular (rough) contact geometries, not only demanding computational resources must be employed for their adequate discretisation, but also convergence issues might arise in standard nodetosegment and segmenttosegment contact strategies, due to corner cases experienced by the contact search algorithms to detect the points in contact. Therefore, regularisation techniques can be employed HPMR ; wriggers , but with the drawback of artificially smoothing out real highfrequency profile features relevant for the physics of contact.
In this work, a relevant extension of the interface finite element with eMbedded Profile for Joint Roughness (MPJR) interface finite element introduced in MPJR is presented to model the fully coupled frictional contact problem between a deformable body and a rigid indenter, with a significant simplification of the contacting interface discretisation, obtained by considering an equivalent smooth interface, instead of explicitly modeling its geometry. The main idea striving the current developments concerns recasting the original geometry of the contacting profiles via a macroscopically smooth interface, which enables the generation of a straightforward meshing with linear finite elements, albeit preserving the actual geometry of the rough profile. This information is stored in terms of its analytical expression through the correction to the initial gap function with the exploitation of the assumption of a rigid indenting profile. The mathematical formulation is detailed in Sec. 2. The advantage of the present approach over alternative meshing methodologies relies on the fact that the shape of the profile, including waviness or even roughness, does not need to be explicitly modeled by the finite element discretisation. In the proposed scheme, the exact interface geometry is embedded into a nominally flat interface finite element through an analytical and exact correction of the normal gap function. In such a way the finite element discretisation and boundary geometry can be taken as nominally smooth, with a significant advantage in terms of regularity of the boundary and simplified finite element discretisation, albeit preserving the exact profile shape in the computations. The full derivation of the proposed finite element is presented in Sec. 3, and has been implemented as a user finite element in the research finite element analysis program FEAP FEAP . The forthcoming contents of the manuscript is arranged as follows. The validation of the proposed model is carried out in Sec. 4, in reference to a benchmark problem, namely an analytical solution of a Hertzian contact problem obtained by assuming semicoupling between the normal and the tangential directions, under different load scenarios. In Sec. 5, the method is exploited to solve cases for which no analytical solutions are available, and that are particularly challenging according to standard finite element techniques. In particular, coupled normal and tangential frictional contact problems with indenters having Weierstrass profiles as boundaries and increasing number of length scales, representative of multiscale waviness, are investigated.
2 Variational formulation for frictional contact problems with embedded roughness
In the present section, the variational formulation which governs the normal and tangential contact of two bodies across a rough interface is presented. First of all, the strong differential formulation is recalled, with the equations describing the mechanics of the two bodies together, along with the HertzSignoriniMoreau conditions for normal contact and the Coulomb law for tangential contact at the interface. Subsequently, the weak form is derived, and the contact conditions are treated. A penalty approach is employed for the normal contact, while a regularised friction law is employed for tangential interactions, in line with finite element procedures wriggers_book . The weak form provides the starting point for the derivation of the interface finite element, which is then presented in the following subsection. The current formulation is valid for D domains, although the proposed framework might be extended to D in a straightforward way.
2.1 General framework
Let us assume that two deformable bodies define the domains in the undeformed configuration defined by the Cartesian reference system . The boundary of the domain can be split into three distinct parts, Fig. 1:

a region where displacements are imposed, i.e. the Dirichlet boundary ;

a region where tractions are imposed, i.e. the Neumann boundary ;

an interface , common to the two bodies, where contact might take place, for which specific boundary conditions must be specified in order to model the stress and deformation fields generated by normal and tangential contact.
As customary, a displacement field , which maps the displacements of the points of from the reference to the current configuration is assumed, identifying and the displacements along and
directions, respectively, and recalling continuous and differentiable functions of the position vector
and time. Further, a small deformation strain tensor is defined as the symmetric part of the deformation gradient,
(with standing for the symmetric part of the gradient operator), which, in standard Voigt notation, reads .At the interface , the configuration of the system is described by the relative displacement field, common for the two bodies, called gap field across the interface, defined as , which is the projection of the relative displacements of both bodies along the tangential and normal directions of the interface defined by the corresponding unit vectors and , respectively. In components, the gap field vector reads .
2.2 Definition of the equivalent contacting geometry
The innovation of the proposed approach lies in the definition of the normal gap . In this framework, regardless of the actual topology, e.g. the one shown in Fig. 2, the contact interface is assumed to be nominally smooth, with uniquely defined normal and tangential unit vectors, while its exact variation from planarity is analytically taken into account as a correction of the normal gap function.
First of all, a local coordinate system is introduced for both the boundaries of the two bodies, being a curvilinear coordinate defining a pointwise correspondence with the coordinates of the same point in the global reference system. Tangential and normal unit vectors and can also be uniquely defined, the latter pointing outwards . If for every point of the contacting interface the indenter’s profile local curvature radius is negligible with respect to the one of the deformable body, i.e. , a smoother line can be set, parallel to and passing through the lowest point of elevation of , Fig. 2. A roughness function is then used to account for the actual profile elevation from the reference datum set in correspondence of . Therefore, the actual profile elevations are described, in the local curvilinear reference system, by the elevation function . The smoothing line is parallel to , and a set of unit vectors and can now be defined, equal and opposite to their counterparts characterised by index .
With this transformation, a zero thickness interface is introduced. It is defined by two distinct lines which at the initial time are perfectly overlapping, and both defined by the function . A relative displacement field can be defined based on the smooth interface of Fig. 3, given by the relative displacements of matching points projected along the normal and tangential direction, e.g. points and of the same figure. The normal gap function defined in this way does not account for the actual complex geometry of the indenter. Its original morphology is stored in the function , which can be combined to with the aim of reconstructing the actual separation between the interfaces, such that it reads:
(1) 
This operation can be performed thanks to the hypothesis of assuming the second body’s material to be rigid, thus guaranteeing that its geometry is kept undeformed and can only undergo rigid body motion. If both the bodies are deformable, the error induced by the present approximation relies in two aspects: (i) the elastic contribution at the interface is locally affected by the geometry regularisation, which has however an average negligible effect for a nominally flat surface where valleys compensate peaks; (ii) the pointwise change in due to local deformation is not accounted for unless an update in the geometry along with the deformation process is considered. Therefore the proposed methodology is expected to provide relatively good engineering approximations to the solution of the contact problem even in the most complex case where both materials are deformable. The real portion of the boundary in contact is now determined by the value of the modified function , in accordance with the condition that contact occurs where . For example, this condition is recovered at assuming that the two profiles make contact only in correspondence of the highest asperity (point of Fig. 3). At the present stage, a negative value for is not admissible since it would imply material compenetration.
2.3 Composite topography
If the contacting bodies can be approximated as halfplanes, e.g if the contact area is small compared to the bulks, the same formulation, i.e. rigid indenter over deformable body, can be interpreted as the solution of a contact problem involving two dissimilar elastic bodies, both potentially characterised by a geometrically complex interface. This problem can be recast in the already proposed framework by defining for the deformable body the composite elastic parameters and , respectively given by:
(2) 
in which and are the original Young’s moduli and Poisson’s ratios, and is the second Dundurs’ constant. A composite topography can be defined ZBP04 ; ZBP07 , combining the geometries of the two contacting boundaries. Again two sets of curvilinear coordinates can be defined, together with normal and tangential unit vectors and , the latter pointing outwards . For each of the two profiles, a smoother line is set, parallel to the average height of the original distribution and passing through its lowest point of elevation. Finally, the actual elevation of is flattened out, so that we have , while the elevation of the indenter is defined as , where now , and are the elevations measured from the respective datum, as shown in Fig. 4. The common coordinate , which is coincident with is here introduced, thus recovering the original case of Fig. 2.
2.4 Governing equations and strong form
The linear momentum balance equation for both and , along with Dirichlet and Neumann boundary conditions on and , can now be recalled for obtaining the strong form of equilibrium for the contacting bodies, enhanced by the conditions for normal contact on , and the classic Coulomb friction law for tangential contact, which determines its further partition as , where the subscripts st and sl stand, respectively for the stick zone, where no tangential relative displacements between the two bodies occur, and the slip zone, where irreversible tangential displacements, i.e. gross slip or sliding, take place wriggers_book :
(3a)  
(3b)  
(3c)  
(3d)  
(3e)  
(3f) 
where denotes the imposed displacement and the traction vector.
2.5 Interface constitutive response
Since not only stick occurs at the interface level, a single penalty approach cannot be used for both normal and tangential responses, thus requiring a distinction for the two different directions. In the normal one, a standard penalty approach has been used, which is characterised by a penalty parameter chosen high enough to reduce compenetration at the minimum. In the tangential direction, a regularisation of the Coulomb’s law is introduced, as in wriggers_book . Given that, the two components of the interface traction vector are:
(4a)  
(4b) 
Equation (4b) is one of the possible regularisations of Coulomb’s friction law which enables a smooth transition from stick to slip condition. Correspondingly, if, on the one hand, it does not lead to a sharp transition from stick to slip domains, it has the great advantage of resolving the continuity issue of the Coulomb’s law at the onset of slip. Finally, the presence of the regularisation variable allows tuning the curve in order to reproduce as closely as possible the real trend of the classic law. It must be remarked that even with a regularised law, the introduction of friction at the interface level is still a source of non linearity for the problem, since tractions do depend on the displacements and therefore on the corresponding solution of the problem.
2.6 Weak form
The normal and tangential contact conditions on , which modify the strong form expressed by Eq. (3) with respect to the classic elastostatic set of equations, also determine the modification of the variational equality representing the weak form of the problem, resulting into the following variational inequality:
(5) 
which is derived from the application of the principle of virtual work together with the application of the contact constraints. In Eq. (5), is the test function (virtual displacement field) for bodies and , which fulfills the condition on . Since body is rigid, its contribution to the integral could be omitted, however it has been taken into account for considering also the more general case of two dissimilar elastic bodies, easily recovering the limiting condition of body being rigid simply setting, in the practical application, a suitable ratio of their respective Young’s moduli. If the contact status is known, e.g if the contact problem is conformal or an active set strategy has been implemented for identifying the contact domain, then Eq. (5) can be recast in the form of an equality by adding the terms related to the normal and tangential contact constraints:
(6) 
where, given the interface constitutive response chosen, the contributions related to normal contact, , and to friction, , read:
(7)  
The displacement field solution of the weak form Eq.(6) is such that it corresponds to the minimum of the energy for any choice of the test functions .
3 Interface finite element with eMbedded Profile for Joint Roughness (MPJR interface finite element) for frictional contact
The numerical solution of the variational problem described by Eq. (6) in the framework of the finite element method requires the geometrical approximation of the two bulks, , and of the interface, , and their discretisation into finite elements:
(8a)  
(8b) 
The bulk has been modeled using standard linear quadrilateral isoparametric finite elements FEAP , even though there is no restriction on the finite element topology, provided that it is consistent with that of the MPJR interface finite element used.
If the contact deformation is small, a conforming discretisation of the bulk and the interface is enough to guarantee the presence of pairs of nodes which are expected to come into contact. This assumption holds also in case of friction and relative tangential displacements, since slip will be infinitesimal too. Given that, a four nodes interface finite element can be introduced, as a special case of a collapsed nodes quadrilateral element. Its kinematics is borrowed from the formulation common in nonlinear fracture mechanics for cohesive crack growth, see OP99 ; PW11b ; PW12 ; Reinoso2014 ; Paggi2015 and then specialised in the present case of frictional contact for the presence of complex surfaces.
The interface element is defined by nodes, each pair belonging to the boundary of one of the two bodies: and to and and to , see Fig. 7. The contribution of the interface to the weak form is expressed by Eq. (7), whose geometrical approximation, together with the finite element discretisation, reads:
(9a)  
(9b) 
In the equations above, the symbol denotes, as customary, the difference between the exact displacement and its finite element approximation, while the subscript spans through the total number of elements employed for the interface discretisation.
The normal gap and the tangential gap are stored in the local gap vector . To evaluate them at every point of the interface element, the nodal displacement vector is introduced, defined as , being and the respective horizontal and vertical displacements of the four nodes. A linear matrix operator is deputed to evaluate the relative displacement between nodes and , and and
respectively. A linear interpolation is then performed by the matrix multiplication with
, which collects the shape functions and . The final step consists in the multiplication by the rotation matrix , for moving from the global to the local reference system, centred in, and aligned to, the interface element. In formulae, it results:(10) 
and the matrix operators take the form:
being , and , the components of the unit vectors normal and perpendicular to the local reference system and . The final tangential gap directly corresponds to , the final normal gap is given by , a correction which accounts for the actual complex surface of the rigid indenter. It has to be remarked that the normal gap is evaluated at every Gauss point, so two times for every interface element. In the discretisation of the contact zone, the number of interface element employed should be high enough in order to adequately sample and reproduce the elevation profile of the embedded contacting shapes. As a rule of thumb, in the case that low order elements are employed, at least elements should be employed for the accurate modeling of each asperity yastrebov:2011 .
At this point it is possible to evaluate tractions and from Eqs. (4a) and (4b), respectively. The contribution of a single interface element to Eq. (6) can now be evaluated. Recalling the traction vector , it is possible to condense Eq. (7), which for a single interface element read:
(11) 
where the variation of the local gap is given, in matrix notation, by:
(12) 
finally, the variation can be set to vanish, and the residual vector can be obtained:
(13) 
which gives:
(14) 
A NewtonCotes integration formula is exploited for evaluating the integral in , which requires the sampling of the integrand at points which coincide to the abscissae of nodes and , obtaining:
(15) 
being the weights and the determinant of the transformation mapping the change of coordinate from the global to the natural reference system. To light up the notation, hereinafter the vector of nodal displacement will be replaced by .
Given the constitutive laws employed, Eq. (14) represents a set of nonlinear transient equations. A full Newton Raphson iterativeincremental scheme is used to solve the system resulting from the discretisation of the rate problem relative to the bulk and the interface. The problem to be solved is therefore the following:
(16) 
The residual vector is linearised introducing the tangent matrix, and the resulting set of linear equations reads:
(17) 
where the superscript denotes the iteration inside a Newton Raphson loop and is the tangent matrix defined as:
(18) 
being the stiffness matrix, the damping matrix and and scalar coefficients that involve the time step and the parameters of the selected time integration scheme. For every cycle, the solution is updated:
(19) 
until the convergence criterion is met. The stiffness and damping matrix resulting from the linearisation of the residual vector for a given iteration read:
(20a)  
(20b) 
where and are the linearised interface constitutive matrices:
(21a) 
and are derived by analytical differentiation if contact is detected, while on the other hand, every term of the matrices is set equal to zero if a positive , i.e. an open gap, is detected.
The application of the proposed approach offers advantages compared to standard FEM discretisation techniques. In particular, there is no need to discretise the interface geometry, since the deviation from planarity can be stored nodalwise with a onetoone correspondence between the original profile and the equivalent interface. This has a twofold advantage: (i) convergence problems of contact search algorithms in case of spiky profiles can be avoided; (ii) the number of finite elements near the boundary can be significantly reduced.
As an illustrative example, let us consider the problem of an elastic block with a wavy boundary, pressed against a rigid flat plane, see Fig LABEL:fig:gmsha. The boundary has been modeled as a Weierstrass profile made of two sinusoidal terms (), see Sec. 5 for details. The equivalent FEM discretisation resulting from the application of the MPJR approach is shown in Fig. LABEL:fig:WM_equiv. In both cases, the same finite element mesh grading as been used, with a characteristic fine mesh size in correspondence of the rough boundary, and a coarse mesh size on the opposite side, being
the overall height of the block. The finite element meshes have been generated using the open source mesh generator GMSH
geuzaine:2009 .The use of the proposed approach is suitable to reduce the number of FEM nodes. At this regard, the reader is referred to the comparison presented in Tab. 1, where rougher profiles with an increasing number of sinusoidal terms have been considered. The gain is due to the fact that the MPJR discretisation does not need to exactly follow the wavy or rough profile as for the standard FEM discretisation, since it embeds the information about the rough topology as a correction to the gap function.

4 Model validation
4.1 Frictional Hertzian contact problem between a cylinder and a halfplane under a monotonic normal load
A semicoupled^{1}^{1}1Hereinafter, the label semicoupled will be referred to a system in which the effect of the tangential tractions over the vertical pressure is neglected, but not the opposite. This hypothesis has been introduced by Goodman goodman1962 . Hertz contact problem is used for validating the model, Fig. 13. This benchmark is found to be particularly suitable for testing the validity of the proposed implementation for two specific reasons. The first underlying motivation, already thoroughly explained in MPJR , is that in spite of its simplicity, the solution of this problem via a standard FEM discretisation of the curved geometry requires very refined meshes, specially at the edges of the contact strips. The second reason is that the present problem yields to a closed solution, thus providing an easy to implement and fast way of comparison and verification. Under the assumptions of purely normal monotonic load and neglecting the influence of tangential tractions on the normal contact traction distribution , the solution can be written as nowell :
(22) 
where is the maximum value of the vertical tractions, and are the extension of the contact radius and slip zone respectively, , and and are the first complete and incomplete elliptic integral of the first kind, respectively. The extent of the slip zone can be evaluated using Spence’s relation spence , which is valid for every initial contact gap defined by a power law relation, and is given by:
(23) 
where is the second Dundurs’ constant. This parameter governs the level of coupling of the system: for , the problem is uncoupled, while its maximum admissible value is . In the limiting case of a contact problem involving a rigid parabolic profile indenting an elastic halfplane, we have , being the Poisson’s ratio of the halfplane. A value of is herein chosen, which corresponds to . While this value is kept constant, four different values for the friction coefficient have been used to investigate the accuracy of the model predictions. The comparison has been carried out in terms of normal and tangential tractions exploiting Eq. (22), comparing our fully coupled numerical approach with the semicoupled analytical solution.
4.2 FEM implementation and results
In the present framework, the actual shape of the contacting profile is embedded in the MPJR interface finite element. The profile elevation is given by its analytical form and evaluated at every NewtonCotes integration point. Since one of the contacting profile, i.e. , is flat, the topography exposed in Sec. 2 simply reduces to the one of : and we have a normal gap in the initial undeformed condition which is given by , being the radius of the cylinder. We analyse the problem under plane strain assumptions. The mesh is structured based on three different levels:

the lower models , which in the present setting is a halfspace, approximated by a circular sector clamped along the curved side and free to undergo vertical displacements in correspondence of the vertical side. An extension of its radius of has been found to be enough for mimicking the elastic properties of a semiinfinite halfplane, under plane strain assumptions. A Young’s modulus and a Poisson’s ratio has been assigned to the standard quadrilateral linear finite elements employed for the bulk discretisation;

the interface is modeled using a single layer of MPJR elements, discretised using elements; The penalty stiffness has been set to , that can be considered as sufficiently high in order to avoid material interpenetration.

finally, the geometry of the indenting cylinder which represents can be replaced by a regular array of quadrilateral linear finite elements, with an assigned elastic modulus
; Neumann boundary conditions are applied as a uniform distribution of vertical pressure
resulting in a unitary vertical force .
The results of the simulations are shown in Fig. LABEL:fig:low_c for , in terms of dimensionless vertical tractions (red curves) and dimensionless tangential tractions (black curves), being the radius of contact related to the semicoupled case and the vertical force applied. For both the analytical semicoupled and the FEM model, five different values of the coefficient of friction, are applied. In the semicoupled case, the distribution of vertical tractions is coincident regardless the coefficient of friction employed, and is highlighted by circle markers, while when coupling is considered, slight differences in the normal tractions can be observed. For what concerns the tangential tractions, an excellent accordance can be observed for all the values of the coefficient of friction employed. The coupling effect observed in the curve related to the numerical approach is in line with the theory, which predicts a stiffening effect that results into an increase of the maximum value of the vertical tractions, compensated by a slight decrease of the contact radius. The slight deviation between the two sets of that can be observed in Fig. LABEL:fig:low_c can be again ascribed to coupling. In fact, the FEM model predicts smaller values for the ratio with respect to the semicoupled approach and the effect can be quantitatively analysed, with results shown in Fig. LABEL:fig:goodman_curve. In this figure, the black solid curve represents the values of for the semicoupled case, as expressed by Eq. (23), corresponding to the dashed black curves of Fig. LABEL:fig:low_c; the red stars are the values for the coupled case, evaluated using an asymptotic solution provided in nowell and the blue circles the outcomes of the simulation. A very good accordance is found between the slip/contact strip width ratio calculated in the fully coupled case and the one obtained by the numerical simulation; the deviation between the tangential tractions of Fig. LABEL:fig:low_c can be justified in this way. The importance of the outcome of Fig. LABEL:fig:low_c lies in the fact that, for some specific instances, care should be taken in neglecting coupling effects, since this could lead to underestimating the magnitude of tractions. As a final remark, it can be noticed that even if a regularised friction law has been used, an adequate choice of the stiffness parameter guarantees results that are very close to the solution based on the Coulomb friction law, which has been exploited by the semianalytical model.
4.3 Hertzian contact problem between a cylinder and a halfplane under constant normal loading and cyclic tangential loading
As compared to the previous test problem, now a downward displacement is imposed on the cylinder, starting from zero and linearly increasing up to a maximum value of . At this point, the vertical displacement is held constant and a tangential load is applied, in terms of a horizontal far field displacement which harmonically oscillates with amplitude . A normalised load history plot is shown in Fig. LABEL:fig:loada, together with the corresponding values of total normal load and tangential load , evaluated as the resultant of the interface tractions. The imposed displacements are also plotted in the load space , in which the black solid curve represents the variation of the tangential load with respect to the normal one, while the blue dashed curves represent the limit of gross sliding. In the present case the load path consists in a curve which is a straight line lying on the horizontal axis, from the origin to point and then becomes a collapsed ellipse with the major axis passing through the points and , Fig. LABEL:fig:loadb. The results in term of traction distributions are shown in Fig. 24. In Fig. LABEL:fig:harm1a, the normal load is linearly increased from zero to its maximum value, and a selfsimilar symmetric central stick area encompassed by two regions of forward slip (positive tangential tractions) and backward slip (negative tangential tractions) develops. Then, Fig. LABEL:fig:harm1b, the tangential load is applied (point of Fig. LABEL:fig:loada and LABEL:fig:loadb), which results in:

an increase in the tangential tractions at the leading edge;

a reduction of the tangential tractions at the trailing edge;

an instantaneous violation of the Coulomb law, which leads to a sudden change from backward slip to stick at the trailing edge, from which tractions start increasing again, but with opposite sign. At the same time, the stick zone shrinks, reaching its minimum at the point of the loading history.
If now the load is reversed, an instantaneous stick zone is created again, which shrinks to a minimum value in correspondence to point , Fig. LABEL:fig:harm1c. The system reaches a stationary state condition after point , but with a steady state value of the stick area which is, for a positive load, sensibly different from the one related to the first application of the load, Fig. LABEL:fig:harm1d. In particular, it retains almost the same extension, but it is shifted towards the center of the contact zone. As a final remark, the system maintains a significant difference between the positive and negative stick areas, even though the steady state is reached after one cycle of loading. It can be noticed that this is directly related to the level of coupling. The current results are obtained with a high level of coupling, corresponding to . The same load history, applied to a system with shows greater similarity between the positive and negative steady stick area, Fig. 29, showing again that the degree of coupling has an important effect on the contact response.
5 Contact between a rough profile and an elastic layer
As shown in MPJR , the strength of the formulation is particularly evident in the case of contact between a rough profile and an elastic layer. In spite of the fact that any kind of elevation field might be taken into consideration in the computational framework, without any restriction, even numerically generated, a Weierstrass profile is herein used as a possible example. The height field is generated by Eq. (24):
(24) 
where in practical applications the summation is carried on up to a certain , thus obtaining a prefractal profile ciavarella ; CIAVARELLA20041247 which consists of the superposition of sinusoidal functions, each of them presenting a decreasing wavelength and amplitude , where and are parameters chosen such that and .
In this section, three different indentation problems are solved, in which the contacting profiles exhibit such heights distribution. Each of them is tested against a rectangular elastic block with a heighttowidth ratio . Such block presents the same elastic parameters employed for the model validation of Sec. 4, i.e. and . Each indenter profile can be considered rigid, and is made of the superposition of one, two, or three terms respectively, according to Eq. (24), shown in Fig. 30. As in the previous section, the lower boundary, i.e. , is flat, and the elevation field reduces to the one of : , with an initial normal gap in the undeformed condition given by . To obtain an adequate sample of the profile heights field, decreasing characteristics interface mesh sizes have been employed, in accordance with the profiles’ shortest wavelength employed from time to time, given as . Following the guidelines of Sec. 3 results in a number of interface finite elements to be employed, depending on the number of terms of the series of , of .
Finally, in order to simulate a contact problem which is indefinite in the direction, periodic boundary conditions have been applied to both the vertical edges of the mesh, at a distance of . A classical ironingtype load history is applied for solving the contact problem. First, a purely normal farfield displacement is imposed, starting from zero up to . Then, a horizontal tangential displacement is applied, linearly varying from zero to times the maximum value of , which is the value that guarantees an incipient gross slip for the single harmonics profile. A full parametric study of the problem should involve a thorough evaluation of the sensitivity of the system with respect to the main governing physical parameters, such as , , , , , and , but this is left for further investigation, since the main purpose is to show the feasibility of treating complex interface problems within the present finite element framework.
5.1 Single harmonics profile in full contact
The case of is shown in Fig. 33. The main difference with the results obtained for the previous test problems is that now we are dealing with an infinitely long profile which makes contact at an infinite set of spots. Under purely normal loading, the vertical tractions present two axes of symmetry, highlighted by dashdotted lines in the figure, which corresponds to axes of antisymmetry for the shearing tractions . Since in this condition shearing tractions must be strictly null in each of these points, they result in having a lower intensity, and the absence of the characteristic backward and forward slip zones which are typical of the Hertzian problem.
Tangential tractions grow in intensity and extension until the full contact condition is reached, corresponding to a value of which is highlighted by the blue dashed curve in Fig. LABEL:fig:WM0a. After that point, since a full stick condition holds, remains constant until the maximum value of the vertical farfield displacement is reached. After that point, the horizontal far field displacement is applied, and the horizontal tractions grow until a condition of partial slip is reached, again blue dashed line in Fig. LABEL:fig:WM0b. As expected, the last point of the interface coming into contact is also the first one which undergoes partial slip. After this point, the state of the system is such that there is an alternation of shrinking stick islands bordered by increasing zones of full slip. When the transient regime ends, a perfect overlapping between and is observed.
5.2 Multiple harmonics profile contact
The addition of a lengthscale in the Weierstrass function has the immediate effect of increasing the peak values of the normal tractions, which are localised in correspondence of the local maxima of the profile, and of reducing the contact area for a given level of the external load: the full contact condition is still achieved, but after a higher number of time steps, see Fig. LABEL:fig:WM1a. The same considerations on the distribution of the tangential tractions can be made also in this case, with the difference that now the contact domain is no longer compact.
In this case, an approximated study of and of the extension of the stick and slip areas can still be possible, exploiting, e.g the CiavarellaJäger theorem, but under the limiting assumptions of uncoupling and equality of the Green functions in and directions. An interesting feature that could be appreciated thanks to having taken into account coupling effects is that, still under purely normal loading, if two asperities, each of them generating a separate contact island, are characterised by a severe gradient in terms of vertical tractions, the less pressed one might experience a horizontal traction distribution which is very far from the one typical of an isolated asperity, i.e. antisymmetric. If the vertical tractions in the leading asperity are high enough, the horizontal displacements generated by them might be so high as compared to the ones generated by the secondary asperity that the latter are negligible, thus resulting in horizontal tractions which are all negative or positive valued from the beginning. With suitable boundary conditions or values of , there might also be a condition of gross slip from the beginning of the contact process.
At the same time, when the second asperity comes into contact, it exerts a stiffening effect on the bulk, which reflects in a decrease of the magnitude of the increment rate of the horizontal displacements towards the high asperity, which in the final place determines a relaxing of the tangential tractions at the level of the leading asperity. This characteristics is depicted in Fig. LABEL:fig:close1 and LABEL:fig:close2
. The tangential tractions over the leading asperity, green dashdotted curves, increase in extension and magnitude as long as the second asperity comes into contact, where they reach their maximum value, blue dashed line. After that moment, they continue growing in extension, since the contact area is still increasing, but they decrease in magnitude, due to the interaction between the different contact islands. Finally, when also the tangential far field displacement is applied, they start growing in magnitude again, gross slip starts, and the transition between full stick and full slip takes place, see Fig.
LABEL:fig:WM1b.The same trend can be observed also for the profile characterised by , Fig. 42, where the same comments apply as well for the increase in vertical tractions, the reduction of the contact area and the evolution of the stick and slip zones.
6 Discussion and conclusion
The proposed formulation provides a way to overcome the shortcomings and the difficulties which are encountered during the solution of contact problem between rough or generically complex profiles, when friction is considered. A standard approach of explicitly modeling the geometry of the interface using higher order interpolation schemes such as Bezier curves, adaptive mesh refinement, or NURBS, is more suitable when the interface consists of regular and smooth profiles, while it can be difficult to exploit or it could be very expensive from a computational point of view when rough profiles have to be analysed.
In this work, the comprehensive and challenging extension to frictional contact problem of the framework which has been set up in MPJR
(denominated as eMbedded Profile for Joint Roughness (MPJR interface finite element)) has been proposed. The fundamental idea was to recast the original geometry of the contacting profiles, obtaining a macroscopically smooth interface which allows for a straightforward meshing with linear finite elements, while the actual geometry is stored in terms of its analytical expression and passed to the system as a correction to the initial gap function, thanks to the assumption of having a rigid indenting profile. The major advantages obtained by exploiting this approach are the use of a low order finite element interpolation scheme, with a significant reduction of nodal degrees of freedom and save of computational time.
The classical benchmark tests which are usually employed for testing the capability of higher order interpolation schemes have been used for validating the model, with excellent results obtained even for relatively coarse interface discretisations. Finally, the method has been successfully tested in relation to more complex scenarios of contact problems involving a Weierstrass profile with multiple harmonics, resulting in a useful tool for the investigation of the behavior of idealised 2D fractal rough surfaces under the non trivial assumption of full coupling between normal and tangential traction fields.
The natural development of the presented interface element, which by itself stems from MPJR , includes its extension to three dimensions and possibly the inclusion of other interface phenomena which could be much more difficult to be analyzed using standard FEM or BEM approaches, such as the interplay of friction and adhesion, or friction and plasticity. The presented implementation also allows for the possibility of more complex friction laws to be used, as, for example, the ones employed in REZAKHANI2020103967 .
Results have highlighted the important role of coupling between normal and tangential contact problems, with a special focus on rough surfaces, which is an open research topic also for precision engineering applications.
Finally, forthcoming studies would encompass the corresponding formulation to geometrically nonlinear effects and prospective coupling of contactinduced fracture events. Such developments are beyond the scope of the present investigation, deserving a careful attention.
Acknowledgements
Authors would like to acknowledge funding from the MIURDAAD Joint Mobility Program 2017 to the project "Multiscale modeling of friction for large scale engineering problems". The project has been granted by the Italian Ministry of Education, University and Research (MIUR) and by the Deutscher Akademischer Austausch Dienst (DAAD). The corresponding author would also like to acknowledge the contribution of the Erasmus+ European project that supported his mobility to Bundeswehr Univerität München (Germany), where he focused on the study of frictional contact problems.
MP would like to acknowledge the support from the Italian Ministry of Education, University and Research (MIUR) to the Research Project of Relevant National Interest (PRIN 2017) XFASTSIMS: Extrafast and accurate simulation of complex structural systems.
JR is thankful to the Consejería de Economía y Conocimiento of the Junta de Andalucía (Spain) contract US1265577Programa Operativo FEDER Andalucía 20142020 and the Ministerio de Ciencia, Innovación y Universidades (Spain) contract PID2019109723GBI00.
Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
References
 (1) J. R. Barber, Contact Problems Involving Friction, Springer International Publishing, Cham, 2020, pp. 41–93. doi:10.1007/9783030203771_2.
 (2) A. Vakis, V. Yastrebov, J. Scheibert, L. Nicola, D. Dini, C. Minfray, A. Almqvist, M. Paggi, S. Lee, G. Limbert, J. Molinari, G. Anciaux, S. Echeverri Restrepo, A. Papangelo, A. Cammarata, P. Nicolini, R. Aghababaei, C. Putignano, S. Stupkiewicz, J. Lengiewicz, G. Costagliola, F. Bosia, R. Guarino, N. Pugno, G. Carbone, M. Mueser, M. Ciavarella, Modeling and simulation in tribology across scales: an overview, Tribology International (2018). (2018). doi:10.1016/j.triboint.2018.02.005.
 (3) A. Klarbring, Examples of nonuniqueness and nonexistence of solutions to quasistatic contact problems with friction, IngenieurArchiv 60 (8) (1990) 529–541. doi:10.1007/BF00541909.
 (4) M. Raous, Art of modeling in contact mechanics, CISM International Centre for Mechanical Sciences, Courses and Lectures 570 (2017) 203–276. doi:10.1007/9783319402567_4.
 (5) M. Cocu, E. Pratt, M. Raous, Formulation and approximation of quasistatic frictional contact, International Journal of Engineering Science 34 (7) (1996) 783–798. doi:10.1016/00207225(95)001212.
 (6) M. Raous, P. Chabrand, F. Lebon, Numerical methods for frictional contact problems and applications, Journal de mecanique theorique et appliquee 7 (1) (1988) 111–128.
 (7) J. Barber, Contact Mechanics, Springer International Publishing, 2018.
 (8) C. Cattaneo, Sul contatto di due corpi elastici: distribuzione locale degli sforzi, Rendiconti dell’Accademia Nazionale dei Lincei 36 (1938) 342–348, 431–436, 474–478.
 (9) R. D. Mindlin, Compliance of elastic bodies in contact, 1949.
 (10) J. Jäger, A New Principle in Contact Mechanics, Journal of Tribology 120 (4) (1998) 677–684. doi:10.1115/1.2833765.
 (11) J. Jäger, Stepwise Loading of HalfSpaces in Elliptical Contact, Journal of Applied Mechanics 63 (3) (1996) 766–773. doi:10.1115/1.2823361.
 (12) M. Ciavarella, The generalized cattaneo partial slip plane contact problem. i—theory, International Journal of Solids and Structures 35 (18) (1998) 2349–2362. doi:https://doi.org/10.1016/S00207683(97)001546.
 (13) M. Ciavarella, The generalized cattaneo partial slip plane contact problem. ii—examples, International Journal of Solids and Structures 35 (18) (1998) 2363–2378. doi:https://doi.org/10.1016/S00207683(97)001558.
 (14) L. E. Goodman, Contact Stress Analysis of Normally Loaded Rough Spheres, Journal of Applied Mechanics 29 (3) (1962) 515–522. doi:10.1115/1.3640599.
 (15) D. Nowell, D. Hills, A. Sackfield, Contact of dissimilar elastic cylinders under normal and tangential loading, Journal of the Mechanics and Physics of Solids 36 (1) (1988) 59–75. doi:https://doi.org/10.1016/00225096(88)900208.
 (16) D. A. Spence, The hertz contact problem with finite friction, Journal of Elasticity 5 (0) (1975) 297–319. doi:https://doi.org/10.1007/BF00126993.
 (17) A. Klarbring, Contact, friction, discrete mechanical structures and mathematical programming, in: P. Wriggers, P. Panagiotopoulos (Eds.), New Developments in Contact Problems, Springer Vienna, Vienna, 1999, pp. 55–100. doi:10.1007/9783709124963_2.
 (18) A. Klarbring, G. Björkman, A mathematical programming approach to contact problems with friction and varying contact surface, Computers & Structures 30 (5) (1988) 1185–1198. doi:https://doi.org/10.1016/00457949(88)901629.
 (19) J. J. Kalker, A Minimum Principle for the Law of Dry Friction, With Application to Elastic Cylinders in Rolling Contact—Part 1: Fundamentals—Application to Steady Rolling, Journal of Applied Mechanics 38 (4) (1971) 875–880. doi:10.1115/1.3408969.
 (20) J. J. Kalker, A Minimum Principle for the Law of Dry Friction—Part 2: Application to Nonsteadily Rolling Elastic Cylinders, Journal of Applied Mechanics 38 (4) (1971) 881–887. doi:10.1115/1.3408970.
 (21) Y. J. Ahn, J. Barber, Response of frictional receding contact problems to cyclic loading, International Journal of Mechanical Sciences 50 (10) (2008) 1519–1525. doi:https://doi.org/10.1016/j.ijmecsci.2008.08.003.
 (22) J. Zhao, E. A. Vollebregt, C. W. Oosterlee, A fast nonlinear conjugate gradient based method for 3d concentrated frictional contact problems, Journal of Computational Physics 288 (2015) 86–100. doi:https://doi.org/10.1016/j.jcp.2015.02.016.
 (23) R. Pohrt, Q. Li, Complete boundary element formulation for normal and tangential contact problems, Physical Mesomechanics 17 (2014) 334–340. doi:10.1134/s1029959914040109.
 (24) K. Willner, Fully coupled frictional contact using elastic halfspace theory, Journal of Tribologytransactions of The Asme  J TRIBOL  TRANS ASME 130 (2008) (07 2008). doi:10.1115/1.2913537.
 (25) J. Leroux, B. Fulleringer, D. Nélias, Contact analysis in presence of spherical inhomogeneities within a halfspace, International Journal of Solids and Structures 47 (2010) 3034–3049. doi:https://doi.org/10.1016/j.ijsolstr.2010.07.006.
 (26) V. Rey, G. Anciaux, J.F. Molinari, Normal adhesive contact on rough surfaces: efficient algorithm for fftbased bem resolution, Computational Mechanics 60 (2017) 69–81. doi:10.1007/s0046601713925.
 (27) V. Popov, R. Pohrt, Q. Li, Strength of adhesive contacts: Influence of contact geometry and material gradients, Friction 5 (2017) 308–325. doi:10.1007/s4054401701773.

(28)
Q. Li, I. Argatov, V. Popov, Onset of detachment in adhesive contact of an elastic halfspace and flatended punches with noncircular shape: Analytic estimates and comparison with numeric analysis, Journal of Physics D: Applied Physics
51 (2018) (2018). doi:10.1088/13616463/aab28b.  (29) C. Putignano, G. Carbone, D. Dini, Mechanics of rough contacts in elastic and viscoelastic thin layers, International Journal of Solids and Structures 6970 (2015) 507–517. doi:https://doi.org/10.1016/j.ijsolstr.2015.04.034.
 (30) J. Zhao, E. A. Vollebregt, C. W. Oosterlee, Extending the bem for elastic contact problems beyond the halfspace approach, Mathematical Modelling and Analysis 21 (1) (2016) 119–141. doi:10.3846/13926292.2016.1138418.
 (31) S. Hyun, L. Pei, J.F. Molinari, M. Robbins, Finiteelement analysis of contact between elastic selfaffine surfaces, Phys. Rev. E 70 (2004) 026117. doi:10.1103/PhysRevE.70.026117.
 (32) L. Pei, S. Hyun, J. Molinari, M. Robbins, Finite element modeling of elastoplastic contact between rough surfaces, Journal of the Mechanics and Physics of Solids 53 (2005) 2385–2409. doi:https://doi.org/10.1016/j.jmps.2005.06.008.
 (33) P. Wriggers, J. Reinelt, Multiscale approach for frictional contact of elastomers on rough rigid surfaces, Computer Methods in Applied Mechanics and Engineering 198 (2009) 1996–2008. doi:10.1016/j.cma.2008.12.021.
 (34) M. Paggi, J. Reinoso, A variational approach with embedded roughness for adhesive contact problems, Mechanics of Advanced Materials and Structures (2018) 1–17. doi:10.1080/15376494.2018.1525454.
 (35) O. Zienkiewicz, R. Taylor, The finite element method for solid and structural mechanics, Vol. 2, ButterworthHeinemann, 2000.
 (36) P. Wriggers, Computational Contact Mechanics, SpringerVerlag Berlin Heidelberg, 2006. doi:10.1007/9783540326090.
 (37) G. Zavarise, M. BorriBrunetto, M. Paggi, On the reliability of microscopical contact models, Wear 257 (2004) 229–245. doi:10.1016/j.wear.2003.12.010.
 (38) G. Zavarise, M. BorriBrunetto, M. Paggi, On the resolution dependence of micromechanical contact models, Wear 262 (2004) 42–54. doi:10.1016/j.wear.2006.03.044.
 (39) M. Ortiz, A. Pandolfi, Finite deformation irreversible cohesive elements for threedimensional crackpropagation analysis, International Journal for Numerical Methods in Engineering 44 (1999) 1267–1282. doi:10.1002/(SICI)10970207(19990330)44:93.3.CO;2Z.
 (40) M. Paggi, P. Wriggers, A nonlocal cohesive zone model for finite thickness interfaces–Part II: FE implementation and application to polycrystalline materials, Computational Materials Science 50 (5) (2011) 1634–1643. doi:10.1016/j.commatsci.2010.12.021.
 (41) M. Paggi, P. Wriggers, Stiffness and strength of hierarchical polycrystalline materials with imperfect interfaces, Journal of the Mechanics and Physics of Solids 60 (4) (2012) 557–572. doi:10.1016/j.jmps.2012.01.009.
 (42) J. Reinoso, M. Paggi, A consistent interface element formulation for geometrical and material nonlinearities, Computational Mechanics 54 (2014) 1569–1581. doi:10.1007/s0046601410772.
 (43) M. Paggi, J. Reinoso, An anisotropic large displacement cohesive zone model for fibrillar and crazing interfaces, International Journal of Solids and Structures 69 (2015) 106–120. doi:10.1016/j.ijsolstr.2015.04.042.

(44)
V. A. Yastrebov, J. Durand, H. Proudhon, G. Cailletaud,
Rough surface contact
analysis by means of the Finite Element Method and of a new reduced model,
Comptes Rendus  Mecanique 339 (78) (2011) 473–490.
doi:10.1016/j.crme.2011.05.006.
URL http://dx.doi.org/10.1016/j.crme.2011.05.006 
(45)
C. Geuzaine, J.F. Remacle,
Gmsh: A 3d
finite element mesh generator with builtin pre and postprocessing
facilities, International Journal for Numerical Methods in Engineering
79 (11) (2009) 1309–1331.
doi:https://doi.org/10.1002/nme.2579.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2579  (46) M. Ciavarella, G. Demelio, J. Barber, Y. H. Jang, Linear elastic contact of the weierstrass profile, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 456 (1994) (2000) 387–405. doi:10.1098/rspa.2000.0522.
 (47) M. Ciavarella, G. Murolo, G. Demelio, J. Barber, Elastic contact stiffness and contact resistance for the weierstrass profile, Journal of the Mechanics and Physics of Solids 52 (6) (2004) 1247–1265. doi:https://doi.org/10.1016/j.jmps.2003.12.002.
 (48) R. Rezakhani, F. Barras, M. Brun, J.F. Molinari, Finite element modeling of dynamic frictional rupture with rate and state friction, Journal of the Mechanics and Physics of Solids 141 (2020) 103967. doi:https://doi.org/10.1016/j.jmps.2020.103967.
Comments
There are no comments yet.