Numerical modeling of adhesion and adhesive failure during unidirectional contact between metallic surfaces

In this work, we developed a finite element modeling approach to study adhesion during unidirectional contact between a two-dimensional plane-strain square and a flat slab. The surfaces were metallic or ceramic, and we analyzed different pairs of materials and their adhesion intensity using a FORTRAN subroutine (DLOAD) connected to a commercial finite element code Abaqus, which provided the surface attractive forces based on the Lennard-Jones interatomic potential using Hamaker constants. We considered adhesive loads during both the approach and separation of the surfaces. During the separation step, we modeled the material transfer between surfaces due to adhesion with respect to damage initiation and propagation at the flat slab. The parameters considered in the simulations include normal load, chemical affinity, and system size, and we analyzed different conditions by comparing the interaction forces during approach and withdrawal. This work also presents: (i) a description of the evolution of energy dissipation due to adhesion hysteresis, (ii) the formation–growth–breakage process of the adhesive junctions and the material transfer between surfaces, and (iii) an adhesive wear map based on a proposed novel equation that correlates the material parameters and material loss due to adhesion. The results indicate that the chemical affinity between bodies in contact is more related to adhesion than the applied load. In addition, the ratio between the material strength and elastic modulus seems to be an important factor in reducing adhesive wear.


Introduction
The adhesion between two surfaces is generally associated with friction. For example, in Bowden & Tabor's plastic junction theory [1], adhesively bonded contacting asperities, so-called junctions, can form, grow, and break during relative sliding [2]. These junctions form in the areas of real contact points, due to surface forces of attraction and repulsion between the atoms and molecules of two approaching surfaces. The adhesive portion is mainly attributed to London van der Waals forces, and may be the source of adhesive failure and wear. As such, adhesion can have an adverse influence on the performance and durability of many mechanical systems, and result in negative economic impacts [3].
According to Hamaker [4], van der Waals interactions between a pair of objects can be obtained by the pairwise summation of the energies acting between all the molecules or atoms in one body with all those in the other body [5]. In this way, the pairwise Hamaker approach represented a first approximation of the van der Waals interactions, which was later complemented by other theories. The Lifshitz theory [6] relates to the absorption properties of real materials, and considers 218 Friction 4(3): 217-227 (2016) macroscopic rather than microscopic quantities, in order to better estimate the Hamaker constants.
The friction force is assumed to be the sum of two contributions: adhesion-related junction-breaking forces and abrasion-related plastic deformation. The first can be estimated by two major parameters-the junctionbreaking shear stress and the junction size-which can be obtained by analyzing the elastic continuum adhesive contact, such as the Johnson-Kendall-Roberts (JKR) [7], Derjaguin-Muller-Toporov (DMT) [8], and Maugis-Dugdale (MD) [9] theories, or by finite element adhesive contact models [10,11,12]. These analyses [712] can also predict the pull-off force, contact force, and contact area [13].
In the JKR model, adhesive forces (attractive tensile forces) are considered only within the contact region, while in the DMT model, adhesive forces are considered only outside of the contact region, with the assumption of Hertzian behavior for deformed profiles [11]. The MD model describes attractive forces based on the Dugdale potential and may be considered a general description of the contact, whereas both the JKR and the DMT models are particular cases. Currently, classical contact models such as the Hertz, JKR, DMT and Maugis models are being used in adhesion experiments [14]. However, these classical contact models do not capture the adhesive failure and/or material transfer between surfaces, as the experimental comparisons are typically based only on the loads acting during contact.
In this work, we modeled the dry adhesive contact between nanoscale, continuum, and homogenous bodies, and considered the case where an indenter applies normal loads on a flat surface. We performed elastic plastic finite element analysis to capture the physical phenomena that promote material transfer between surfaces, and evaluated the effect of system size and material properties on the adhesion forces. We estimated crack generation, which leads to material transfer between surfaces, and which provides information regarding the critical conditions that promote wear due to adhesion.

Adhesion forces and Hamaker expressions
The Lennard-Jones potential, whose most common form is given by Eq. (1), is a simple physical model that approximates the interaction between a pair of neutral atoms or molecules. 12 6 ( ) 4 w r r r In Eq. (1), ( ) w r is the Lennard-Jones potential energy function, r is the distance between particles,  is the depth of the potential well, and  is the distance at which the inter-particle potential is zero. In this equation, the first term describes the short-range repulsive interactions and the second, the negative term, refers to the long-range attractive interactions. According to Eq. (1) and Fig. 1(a), when two approaching surfaces or particles are placed closer than within a few nanometers, they will be subjected to attraction forces. When differentiating Eq. (1) with respect to the separation distance, r, we obtain the graph shown in Fig. 1(b), which indicates that attraction will continue until the separation distance is equal to 1/6 2  , when repulsion forces may be initiated. For separation distances below 1/6 2  , contact mechanics can be used to describe the mechanical behavior. In this situation, the pair potential between two atoms or small molecules can be simplified and considered   as purely attractive, for example, considering only the van der Waals forces. In this case, Eq. (1) can be rewritten as: If the interaction between two flat surfaces (one being infinite) is non-retarded and additive, and the separation distance between surfaces is much lower than the object thickness [15], Eq. (2) can be integrated based on the energy of the van der Waals interactions between one molecule of the first surface and all the molecules of the second surface. We do so by integrating the molecular density throughout the solid volume, resulting in an interatomic van der Waals pair potential, per unit area, as represented in Eq. (3) [15]: Differentiation of Eq. (3) provides the adhesion force (per unit area), as follows: where A is the Hamaker constant for the material pair.
In this work, in order to avoid the adhesive pressure becoming infinite, r never goes to zero after considering that 2 where h is the clearance between surfaces and  is the atomic radius of the element present in the surface. For two macroscopic phases, 1 and 2, which interact across a medium 3, A can be estimated based on the optical and electric properties of a material using McLachlan's equation [15]: where the input parameters are the dielectric constants  ( ) and the refractive indices (n), and e  is the absorption frequency. For two dissimilar materials, the Hamaker constant may be estimated in terms of the geometric mean of the Hamaker constant of each material [16], that is: 12 11 22 or as an alternative [16]: 11 22 12 11 22 2A In this work, we obtained Hamaker constants for the same material in both contacting surfaces from the Refs. [15,17] and we used Eq. (6) to compute the Hamaker constants for different materials in contact.
The numerical approach considers only the adhesive part of the potential and we modeled the repulsive part of the interatomic potential using the standard formulation for normal contact, as implemented in the FEM software Abaqus, which is referred to as "hard" contact. For this behavior, a contact constraint is applied when the clearance between the two surfaces becomes zero. In the contact formulation, there is no limit on the magnitude of the contact pressure that can be transmitted between the surfaces. The surfaces separate when the contact pressure between them becomes zero or is negative and the constraint is removed.

Adhesive failure and material transfer simulation using XFEM
Adhesive forces produce stresses and strains that may locally exceed the material strength and lead to material failure and, consequently, lead to material transfer between the surfaces in contact. In this work, we applied the extended finite element method (XFEM) to reproduce the material cracks that arise near the contact region and promote material transfer and wear. With this technique, a crack is nucleated based on a fracture initiation criterion, which we defined as the maximum principal stress [18]. The crack propagates according to a damage evolution criterion based on the energy release rate. In other words, crack initiation is based on the maximum tensile stress value of an element in the system. When the maximum principal stress reaches the predefined tensile strength of the material, a crack is initiated and can propagate. In the case of a homogenous stress field, the crack will propagate in a direction perpendicular to the maximum principal stress. Crack propagation is based on strain energy release, which means that the extension of the crack depends on the balance between the surface energy and strain energy release rates.

Model description
We conducted a finite element analysis to investigate the contact problem of a linear elastic square punch indenting an elasticplastic deformable slab. To this end, we used the finite element solver Abaqus to model the elastic adhesive unidirectional contact of a deformable (square or rectangular) asperity with a plane-strain surface. We coupled the two-dimensional plane-strain analyses with an ad hoc user FORTRAN subroutine designed to calculate the adhesion forces, which we introduced in the system as forces acting on the surfaces as a function of their separation distance. We used static implicit modeling and simulation to determine the stresses, strains, and displacement response and to evaluate the material failure. The system geometry, as presented in Fig. 2, indicates dimensions that are compatible with the contact of a nanometric asperity indenting a flat plane surface. In contrast to some other analytical adhesion models, such as the DMT model [8], in this work, we considered the interaction between the indenter lateral faces and the slab to be negligible.
In the simulations, we considered two different indenter edge lengths, as illustrated in Fig. 2. The model mesh was composed of 1,092 nodes and 1,000 elements in configuration A, and 982 nodes and 900 elements in configuration B.
With respect to the boundary conditions, the slab was fixed at the bottom line, and neither lateral line could move in the x-direction. The indenter was allowed to move only in the y-direction. The initial separation distance between the indenter and the slab was 0.3 nm. We imposed normal load by two penetration depths of the indenter over the slab, 0.1 nm and 0.2 nm, which are denoted as 1 and 2, respectively. Penetration depth is defined as the distance that the bottom edge of indenter moves below the initial position of the slab's upper edge. We applied this movement to the top line of the indenter. We restricted the numerical increment size to always being less than 0.01 nm/increment, in order to avoid abrupt variations in the adhesion forces.
We chose the penetration depth values with respect to the dimension of an atomic diameter. Despite the very shallow penetration depths, we emphasize that the analyzed system has nanometric dimensions, since adhesion forces arise at this length scale.
A molecular dynamics formulation would be more appropriate for small system simulations. However, a continuum approach at this scale is more simply expanded to bigger systems than would atomistic formulations.
To avoid inconsistencies in the physics of the system, since fracture can occur in the elements of a nanometric system and result in the detachment of extremely thin layers, choosing an element size bigger than an atomic diameter is appropriate. However, to avoid theoretical inconsistencies when using a continuum model, the size of the elements in the mesh must be chosen to be fine enough to accurately capture the investigated properties and be independent from any atomistic approach. In addition, while the continuum theory should only be used down to a certain length scale, once a continuum theory is chosen, the mesh size should be as small as possible to insure accurate results irrespective of the atomic distances.
To address this issue, we first performed a mesh convergence analysis, and found that the use of an element size smaller than 0.3 nm allows for convergence in the numerical results. Therefore, we chose an element size of 0.25 nm, which is also larger than an atomic diameter. Figure 3 shows the mesh convergence analysis results. We can observe singularities for element sizes around 0.5 nm, and a convergence in the results for the stress field and the amount of material transferred between surfaces for element sizes smaller than 0.3 nm. In this situation, we can continue to use an element size bigger than an atomic diameter, since it meets the requirements of the continuum approach and also allows for comparison with atomistic systems.
For the indentation depth and size of the contact region in the adhesion contact models, we chose values for which the bulk and surface were comparable with molecular dynamics simulations.
The most important inputs to the finite element model include: (i) the interfacial geometry; (ii) the mechanical properties of sapphire, copper, and H13 tool steel; and (iii) the adhesion parameters (Hamaker constants). Using 24 combinations, we systematically and parametrically varied the model inputs to obtain the expected variability, as described in Table 1.
To describe the mechanical behavior of the indenter and the slab, we used an elasticplastic material model. We obtained the elasticplastic constitutive properties of the materials from the Refs. [19−21], as presented in Table 2. We imposed the condition that only the slab elements could fail and be transferred to the indenter, allowing no material transfer from the indenter to the slab.
In this work, we adopted a single procedure in the XFEM method, although it has been recognized that further improvements are possible in the description of copper and other ductile materials [22,23].
XFEM analysis requires additional material information in the enriched area of the crack. The term enriched area is defined by Abaqus as the region  where a crack can nucleate and propagate [24]. We selected maximum principal stress criteria for the damage initiation criteria. For all evaluated materials, we used tensile strength values (as shown in Table 2) as the limiting maximum principal stress. We used fracture energy for the damage evolution criteria, which we estimated using the following Irwin fracture condition: where G I is the fracture energy, K I is the fracture toughness, and E is the elastic modulus.
The Hamaker constant, A (in Joules), represents the strength of the van der Waals interaction forces between macroscopic bodies, which can be roughly defined as a material property. Table 3 lists the values used in our numerical simulations, which were obtained either from the literature or from Eqs. (6) and (7). We used iron oxide properties to define the adhesion  parameters (Hamaker constant) in the contact of the H13 steel surfaces, since iron and steel surfaces are usually found to be coated with oxide layers. At the same time, we considered the mechanical properties of H13 steel for the material below the surface.

Ad hoc subroutine for adhesion forces calculation
Based on the potential and force expressions in Eqs. (3) and (4), we developed, compiled, and connected a user-defined load FORTRAN subroutine (DLOAD) [24] to the commercial finite element code Abaqus to calculate the adhesion forces during the approach and separation of the surfaces. The subroutine iterates at each time increment and for each element searches for the closer element in the opposite surface. It updates the surface separation at each surface integration point and uses this value to feed Eq. (4), thus providing the surface forces in the elements at each surface. These surface forces are then applied to the edges of surface elements. At this stage, we used the Hamaker constants presented in Table 3.
This approach allowed for the development of a fully coupled strainadhesion forces model, since adhesion forces are computed as a function of the separation distances between surfaces, while taking material deformation into account.
During unloading, the simulation focused on adhesive failure. We modeled the material transfer between surfaces due to adhesion by the damage initiation and propagation criteria in the XFEM, such that we could investigate the mechanical parameters that influence these phenomena. We calculated the variations of the interaction force during approach and withdrawal, and the dependence of pull-in and pull-off force on the surfaces' approach and separation. Figure 4 shows a simplified flowchart describing the iterative analysis process and the coupling between this subroutine and the FEM software.

Results and discussion
In this work, we estimated the wear related to adhesion phenomena in dry contacts based on the amount of material transfer from the slab to the indenter. We then evaluated this transfer for each of 24 different combinations of surface pairs and, based on the obtained results (Table 4), proposed an adhesive wear map.
We also evaluated the morphology of adhesive failure and some energy considerations. The results are organized in subsections below. Figure 5 shows the final configuration of the selected systems. These examples represent the four typical behaviors we observed from our evaluations of the adhesive failure and material transfer for 24 input combinations. Table 5 presents the observed configurations for each of the 24 cases, with the results organized according to the amount of material transferred between the surfaces, i.e., (i) no material transfer and  Our results indicate that the amounts of material transferred from one surface to another depend on the slab elastic modulus and other mechanical properties. In addition, materials that have lower hardness experienced higher material transfer. Increasing values of the Hamaker constant as well as increasing penetration depths also contribute to failure and material detachment.

Types of adhesive failure
The main reason for the detachment of large layers of slab material seems to be increasing penetration depth, which in association with increasing contact point areas, leads to the highest values of material loss.   We note that, for the oxide-covered H13 slab, a simultaneous increase of load (penetration depth), indenter size, and chemical affinity did not produce any material transfer. This means that, in these conditions, higher mechanical property values, such as the ratio of the tensile strength to the elastic modulus, can minimize the effects of adhesion.

Adhesion hysteresis
With respect to the contact loading and unloading cycles, an adhesion hysteresis may occur if there is a difference between the forces during the loading and unloading stages. This phenomenon is not captured by classical contact theories, such as Hertz, JKR and Maugis-Dugdale, and has been attributed to dynamic effects, such as those due to material viscoelasticity and ambient and surface moisture, or to materialspecific chemical reactions [25]. However, many quasistatic contact measurements also display hysteresis and it is possible to have depth-dependent hysteresis during perfect elastic contacts.
The hysteresis observed in this study was generated mainly by plastic deformation and also by the bonding between surfaces. Figure 6 shows an example of hysteresis behavior in the reaction force observed in the numerical results. These results are similar to those captured in molecular dynamics simulations [25] and experimental measurements [26], but were captured here using a numerical approach based on continuum mechanics.
It is possible to identify two different behaviors in the hysteresis curves. In Fig. 6(a), showing simulations that did not result in a material transfer between surfaces, we see an unloading force peak slightly greater than the peak related to the loading stage (case 01, for example). In situations in which material transfer was observed, for example in case 14 ( Fig. 6(b)), the unloading peak force is less apparent or even absent and is always smaller than the loading peak.
One possible explanation for the observed hysteresis phenomenon in this work is that adhesive failure can promote energy loss due to crack propagation in the material, leading to less energy being required to separate a pair of bodies in contact. Therefore, an increase in the energy loss due to material failure is equivalent to an adhesive toughening of the contact interface. This observation seems to be in agreement with the work of Zheng and Zhao [27], in which adhesion hysteresis is built into the assumption that different amounts of work must be done when bringing materials into and out of contact.
Due to the very large volume of results obtained in the numerical simulations conducted in this work, we show only the results of cases 01 and 14 in Figs. 6(a) and 6(b), respectively. However, we note that the results of the other cases are similar when considering the groups presented in Table 5. Therefore, in all cases for which a layer is detached from the surface, we observed a low peak in separation, as shown in Fig. 6(b), and different behavior when adhesive wear was not observed, as in Fig. 6(a).

Adhesive wear and material transfer
Our attempts to organize the results obtained in this work led to our developing a definition of the parameter adhesion  to estimate the intensity of adhesive  wear. This parameter is presented in Eq. (9), where depth h is the penetration depth, L real is the real area of contact, 12 A is the Hamaker constant, E is the elastic modulus, and failure  is the material stress strength. This parameter is analogous to Archard's equation [28], since it relates mechanical properties and system characteristics to evaluate wear, but is dedicated to the analysis of adhesion and, for now, is restricted to the nanoscale.

 
depth real 12 We generated Fig. 7 based on this expression, resulting in our proposed wear intensity map for adhesion. By analyzing Fig. 7, we can identify three different regions. We can relate the first region, for adhesion  values between 0 and 25, to the absence of wear due to adhesion between surfaces. Values between 25 and 125 can be associated with mild adhesive wear, and adhesion  values greater than 125 indicate severe adhesive wear.

Conclusion
From the simulations conducted in this work, we can draw the following conclusions:  Using a numerical approach based on continuum mechanics, we produced results that are consistent with those obtained previously in molecular dynamics simulations [25,29].
 Using this numerical approach, we could reproduce adhesion hysteresis.
 Wear due to adhesion is affected by chemical affinity, as well as by the normal load (or indentation depth), contact area, and, significantly, by the mechanical properties of the materials in the contact pair.
 There seems to be a strong correlation between the amount of material transferred between surfaces and the proposed parameter adhesion  . Based on our results, we identified three distinct adhesion groups, which we illustrated using an adhesive intensity wear map with adhesion  as the abscissa.
Open Access: The articles published in this journal are distributed under the terms of the Creative