An adaptive mesh refinement approach for solving non-Hertzian elastic contact problems

Semi-analytical methods are a common way of solving non-hertzian contact problems when designing mechanical components. These methods require of the discretization of the domain into a set of pressure elements and their accuracy and computational cost are related to the number of elements in which the domain is discretized. But, while the accuracy increases as the pressure element mesh is refined, the computational cost increases quadratically with the number of pressure elements. So in the great majority of the cases, a commitment between accuracy and computational cost must be achieved. In this work, a new approach has been developed to improve the performance of semi-analytical methods for solving contact problems. This approach uses an adaptive mesh refinement strategy, based on the quadtree decomposition of the domain. As a result, the computational cost decreases, while the accuracy of the method remains constant.


Introduction
The contact stress analysis plays an important role during the design process of several mechanical elements like roller bearings, gear transmissions, etc. In order to accomplish this contact analysis, the socalled contact problem must be solved to obtain the following relevant information: (1) The contact area, that involves the determination of the size, shape, and location of the contact areas in each one of the contacting bodies. (2) The contact stresses, that involves the determination of the contact pressure distribution on the surface of the bodies and also the stress distribution underneath the surfaces. (3) The deformation of the bodies produced by the contact pressure.
Different approaches have been used to solve the contact problem when designing mechanical components. Some of them are numerical methods based on the discretization of the domain, such as the finite element method (FEM) [2,4]. Other type of approaches use analytical solutions [3], generally based on Hertz contact theory [10], to determine the stress and the displacement fields produced by the contact. Compared to the FEM approach, it can be said that the analytical approaches are more efficient in terms of computational cost, but they have severe applicability limitations, since the assumptions of the Hertz theory must be met. On the other hand, the FEM approach can overcome these limitations, but at a much higher computational cost. Semi-analytical methods (SAMs) [8,13] can be considered as an intermediate solution to the contact problem. Many of these methods solve the contact problem by means of influence coefficients that relate contact pressures to surface displacements. The influence coefficients may have different nature, but in the great majority of the cases, they are determined using the Boussinesq-Cerruti analytical solution for point loads acting on an elastic half-space [12,14]. This method requires of the discretization of the domain in n pressure elements and then, the total displacement at any point is determined by superposition of the displacements produced by contact pressures acting on each pressure element of the domain. These methods are potentially more efficient and faster than the FEM approach and they allow to overcome some of the limitations of the analytical approaches in nonhertzian contacts.
For these reasons, SAMs have been commonly used in machine design. Harnett [7] and de Mul et al. [16] performed bearing contact calculations using SAMs. Jin [9] and Pascal [18] used these methods to determine contact pressure distribution in wheel/rail contact. Finally, SAMs have also been used in the field of gear transmissions, where relevant works have been carried out by Sheveleva [20], Wu [22] and Guilbault [6].
Similarly to the FEM approach, the computational cost of the SAMs depends on the desired accuracy when determining the contact pressure distribution and the true contact area. Kalker [11] demonstrated that the computational cost of these methods is proportional to n 2 and stated that the accuracy of the SAMs is defined by the number of pressure elements in which the true contact area is discretized, specially in the areas close to the border of the contact. This implies that, in the great majority of the cases, a commitment between accuracy and computational cost must be adopted.
When both shape and location of the true contact area are known in advance, the efficiency of the method can be maximized by discretizing an area similar to the true contact area. But when the true contact area is unknown, it is difficult to optimize the efficiency of the method, since the whole potential contact area must be discretized to consider any possible shape and location of the true contact area. In consequence, there could be many pressure elements in the discretization out of the true contact area, what causes a loose in the efficiency of the method.
These difficulties could be partially overcome using adaptive mesh refinement strategies. These techniques have been previously used to improve the efficiency of the numerical methods based in the discretization of the domain, specially in FEM procedures [17,21]. However, no previous use of adaptive refinement has been found in the literature for the solution of nonhertzian contact problems through the influence coefficient method.
In this work, the influence coefficient method is used to define an algorithm to solve non-hertzian contact problems, which performance is improved by using an adaptive mesh refinement strategy. The proposed algorithm is based on the quadtree decomposition of the domain, and the mesh refinement is performed based on the rate of change of the contact pressure, that is a discrete estimation of the gradient of pressure. A parametric study of the performance of this new approach is performed and illustrated with point and line contacts.

Description of the contact model
In this work, a contact model based on influence coefficients has been used to solve the frictionless elastic contact between two bodies, under the assumption that both of them can be approached to elastic half-spaces in the vicinity of the contact. This approach provides a powerful method to predict the contact pressure distribution that appears when two bodies become in contact. It is based on the superposition of the Boussinesq solution for a point load applied on an elastic half-space, and its main ideas are explained below.
Let us consider two bodies 1 and 2 in its undeformed contact position, contacting at the initial point of contact O L (Fig. 1a). In this point, a common tangent plane P is defined, that is assumed to be so close to the surface of the bodies in the vicinity of the contact area that the deformation of the surfaces of both bodies may be referred to it in the linear small strain theory of elasticity. A Cartesian coordinate system is defined with origin at O L , being the local axis Z L normal to the plane P and pointing inward the body 2.
Consider a generic point Q in the plane P, whose position is defined by the vector rðx L ; y L ; z L Þ, being z L ¼ 0. The gap between the two bodies at point Q, measured along the Z L axis, is denoted by the function BðrÞ.
The two bodies are pressed together in absence of friction by the effect of the force F T (Fig. 1b), causing an approximation between them that is denoted by d. This distance can be defined as the approach, in the direction of the Z L axis, between two reference points (each one associated to each body) that are placed in regions of the bodies that are not affected by the local deformations produced by the contact.
Since penetration is physically inadmissible, a contact pressure distribution pðrÞ is generated in the contact area A, which produces a normal elastic deflection xðrÞ in the bodies. Both magnitudes are related by the Boussinesq Eq. [14]: where r 0 is the position vector of dA, m is the Poisson coefficient, G is the shear modulus and superscripts 1 and 2 refer to each of the bodies involved in the contact.
As stated by Kalker [11], the true contact area and the pressure distribution may be determined minimizing the total complementary energy V (Eq. 2) under the condition that contact pressures are positive everywhere.
To enable the numerical solution, the potential contact area is discretized into a set of n rectangular pressure elements D j (j ¼ 1::n) of area A j , as shown in Fig. 2, where the contact pressure on each element is assumed to be constant (pðrÞ ¼ p j ). The position of the centroid of each element is defined by vector r i . Under this discretized domain, Eq. 1 can be rewritten as: where f j;i is the influence coefficient of pressure element D j over the centroid of element D i . If the pressure elements have a rectangular shape, the magnitude of f j;i may be determined as described by Johnson [10]. Finally, if BðrÞ is assumed to be constant over each pressure element (Bðr i Þ ¼ B i ), Eq. 2 may be rewritten as: For a given value of d, the solution to the contact problem can be found, in terms of the contact pressure distribution, by minimizing V under the following conditions: However, in the great majority of the cases the solution of the contact problem is required for a given value of F T instead of being required for a given value of d. In those cases, an additional force equilibrium constraint can be added to the minimization problem, which is given by: In this way, not only the contact pressure distribution can be calculated for a given value of F T , but also the corresponding approximation d between the contact bodies. The true contact area is defined, within the precision of the mesh size, by the boundary between the zero and non-zero pressures. The size of the true contact area A is calculated as the sum of the areas of the pressure elements A i where p i [ 0.
If any of the contact bodies has finite dimensions, the assumptions of the proposed contact model are no longer fulfilled, because that body cannot be approached to an elastic half-space, and in consequence, an erroneous solution of the contact problem may be obtained. In these cases, the mirroring method proposed by de Mul [15] combined with the overcorrection factor proposed by Guilbault [5] can be used to take into account the influence of the free surfaces over the contact. This method involves the calculation of additional influence coefficients of pressure elements which are mirrored respect to the end planes of the bodies, multiplying the computational cost by ðM þ 1Þ, being M the number of finite dimensions taken into account.
3 Quadtree decomposition of the domain According to Samet [19], the basic concept of the quadtree is to enclose the domain of the problem C into a containing cell, usually a square, which is denoted as the root of the quadtree, as shown in Fig. 3a. This cell is, then, subdivided into four sons of the same size (Fig. 3b), one in each direction: North-West (NW), Nord-East (NE), South-West (SW) and South-East (SE). Each one of these cells is subdivided recursively until a stopping criterion is reached, which may be based upon the local geometry of the domain or in user-defined parameters (Fig. 3c, d). The information related to the quadtree decomposition of the domain is stored in a hierarchical tree structure, as shown in Fig 3e. For every cell, references to its ancestor and sons are stored. This kind of structure ease the performance of several operations, such as the neighbour finding in a defined direction, that will pay an important role in the proposed method.
Each corner of a cell is called vertex. The level of a cell j in the structure is denoted by L j , and represents the number of divisions performed from the root of the quadtree. According to this definition, L j is also related to the relative size of the cell inside the quadtree structure and the degree of mesh refinement that this size represents. Given the size of the root cell of the quadtree, the size of any cell can be determined if its degree of refinement L j is known. The root cell of the quadtree is usually denoted by level 0.
Any cell that is not subdivided anymore is a leaf cell (displayed in grey in Fig. 3e), while subdivided cells are referred to as non-leaf cells.

Contact algorithm with adaptive refinement
The application of the contact model described in Sect. 2 requires the discretization of the potential contact area in a set of n pressure elements D j . As usual in numerical methods based on the discretization of the domain, the election of the number of pressure elements in which the domain is divided involves a commitment between accuracy and computational cost. Kalker [11] stated that the computational cost of the algorithm is proportional to n 2 , that is the number of influence coefficients that need to be calculated to solve the contact problem. He also argued that the accuracy of the model to describe the true contact area depends on the mesh refinement, specially in the areas close to the border of the contact area. Consequently, an improvement of the accuracy of the results necessarily implies an increment of the computational cost.
When the size and position of the true contact area are known in advance, the efficiency of the model can be easily maximized by discretizing the portion of the domain that fits the true contact area. But when they are unknown, the whole potential contact area must be discretized to take into account any possible shape and location of the true contact area. Then, it is usual to use a uniform mesh for the whole domain [5,22], being more or less dense depending on the desired accuracy and on the capabilities of the computer used to solve the contact problem. In these cases, the efficiency of the method gets worse as the true contact area gets smaller compared to the potential contact area. And this effect is accentuated when the contact bodies have finite dimensions, since extra operations are required to take into account the influence of the free surfaces over the contact area, as described in Ref. [15].
To cope with that problems, a contact algorithm with adaptive mesh refinement has been developed. This algorithm uses the quadtree decomposition of the domain (described in Sect. 3) as discretization to solve the contact problem. In this case, the domain to be enclosed by the root cell of the quadtree is the interference area ðCÞ, which is defined as the boolean intersection of the projection of the bodies on the plane P, as shown in Fig. 4. All the leaf cells of the quadtree are considered pressure elements.
To maximize the efficiency of the proposed algorithm, it is important to minimize the number of leaf cells of the quadtree (pressure elements) located outside of the interference area. This can be achieved by ensuring that the interference area is enclosed by a root cell of the quadtree coincident with the minimum bounding rectangle (MBR), defined as the minimum rectangle that contains every point in the region [1]. If the interference area has a predominant direction, multiple quadtrees may be joined together in that direction, to avoid undesired aspect ratios of the resulting pressure elements.
The use of a quadtree data structure offers two interesting features to this algorithm. In first place, the  that is a parameter that describes the size of the elements of the initial uniform pressure element mesh.
The algorithm starts determining the common tangent plane P, where a local Cartesian coordinate system is defined, being the Z L axis normal to the plane P (step A1). The boundaries of the contacting bodies are normally projected onto the plane P to determine the interference area C (step A2).
The interference area is enclosed by the root cell of one (or more) quadtree that is determined using a general purpose algorithm to find the MBR [1]. The local Cartesian coordinate system is rotated so the X L and Y L axes become aligned with the principal directions of the MBR. The root cell of the quadtree is, then, recursively subdivided until the desired initial level of uniform mesh density L ini is reached for all the cells of the quadtree (step A3), obtaining a uniform mesh of pressure elements. All leaf cells of the quadtree are considered pressure elements D j for the initial iteration of the algorithm. All the pressure elements are marked with the flag K i ¼ TRUE indicating that their properties (normal gap B i and influence coefficients f j;i ) are not computed yet.
Then, an iterative process starts whose first step is the determination of the normal gap B i for all the pressure elements (step A4). The influence coefficients f j;i of those elements are also determined (step A5) following the equations described by Johnson [10].
Finally, the solution to the contact problem is found by solving the linear system of equations given by Eqs. 5 and 6 (step A6). As explained by Wu and Tsai [22], since the interference area overestimates the true contact area, negative values of p j will be obtained when solving this system of equations. In order to satisfy the conditions given by Eq. 5, an iterative process starts, in which negative values of p j are set to zero and the corresponding pressure elements D j are discarded for the next iteration. Then the system of equations is solved again, and this iterative process is repeated until the condition p j ! 0 is satisfied. As a result, the contact pressure distribution p j and the corresponding approximation between the contact bodies d are obtained.
The flag K i is defined as FALSE for all the elements present in the discretization, indicating that the properties of these elements have already been computed. Then, the algorithm to determine the elements to split (described in Sect. 4.1) is called (step A7), returning an array with the indexes of the elements that must be split. These elements are split (step A8) and the quadtree data structure is updated with the information of the new elements, which are marked with the flag K i ¼ TRUE, indicating that their properties are not computed yet. If no new elements are created, the iterative process finishes and the contact results are displayed (step A9). In contrast, if new elements are created, the iterative process starts again (step A4), and it is repeated until no new elements are created.
The main advantage of this algorithm is that only the normal gap and the influence coefficients related to the new elements are computed for each iteration, decreasing the global computational cost of the algorithm. The number of influence coefficients calculated by the algorithm ðN f Þ can be determined afterwards using the following equation: where t is the number of iterations performed by the algorithm, n ðiÞ is the number of elements in iteration i, n newðiÞ is the number of new elements in iteration i and n oldðiÞ ¼ n ðiÞ À n newðiÞ .

Algorithm to determine the elements to split
An adaptive mesh refinement may be based upon several criteria, such as the local geometry of the domain or the rate of change of physical magnitudes. In a contact problem under the small strain domain, the contact area is usually small (compared to the dimensions of the contacting bodies) and the gradient of pressures is usually high in part of the contact area. Moreover, for a given level of accuracy, the approximation to the pressure distribution with discrete elements of constant pressure requires small elements in areas where the gradient is high, but the elements can be larger where the gradient is low. For this reason, the gradient of the contact pressure is used as criterion for refinement in this work. But, since the algorithm works with a discretization of the domain, the gradient is estimated as the rate of change of the contact pressure ðu j;i Þ between two adjacent elements D i ; D j . Finally, to make the parameter dimensionless, it is defined as the absolute difference of the contact pressure of element D i and a neighboring element D j divided by the average value of these pressures: The obtained value for u j;i is compared to an arbitrarily defined value u max (representing the maximum allowed rate of change of contact pressure between adjacent elements) to decide if the related elements should be split or not. Therefore, when u j;i [ u max both pressure elements are marked as candidates to be split. However, in the border of the contact area the pressure function is not differentiable. Then, according to Eq. 8, the rate of change between an element D i that is within the contact area (having contact pressure p i [ 0) and an adjacent element D j that is outside of the contact area (having no contact pressure, p j ¼ 0) is u j;i ¼ 2 regardless of the value of p i . As the value specified for u max will usually be lower than 1, it can be concluded that this refinement strategy based on the rate of change of the contact pressures will refine the mesh at the boundary of the contact area endlessly.
In order to limit the number of iterations performed by the algorithm, an additional stopping criteria based on the minimum size allowed for a pressure element is included. As mentioned before, the level L j that a pressure element occupies in the quadtree structure is related to its size, so limiting the former will also limit the latter. This limit is defined by an user defined parameter L max , referred to as the maximum degree of mesh refinement.
It is important to point out that u max is a target value and may not be always reached. If the maximum degree of mesh refinement L max is reached for all pressure elements before the target value for u max is achieved, the algorithm will finish. This guarantees that the maximum degree of mesh refinement will be reached at the border of the true contact area, so the precision in which the border of the contact area is determined is defined beforehand by L max , regardless of the value selected for L ini .
Other criteria to perform and stop the mesh refinement could be implemented, but depending on the involved variables it could lead to some inconveniences, like a higher time of computation or a poorer description of the contact pressure distribution. Since the contact pressure is the primary result of the approach, having a good description of the pressure distribution is essential to get realistic values of the derived variables (like displacements). For this reason, a refinement criterion based on pressure is highly recommended here.
The main routine of the algorithm to determine elements to split is shown in Fig. 5. The following input information is required by the algorithm: (1) An array containing the contact pressures p j associated with every pressure element of the current discretization. (2) The quadtree data structure.
(3) The maximum degree of mesh refinement ðL max Þ. (4) The maximum allowed rate of change of contact pressure between adjacent elements ðu max Þ.
The algorithm starts defining the flag K i as FALSE for all the elements present in the current discretization. The flag K i indicates when a pressure element must be split, so in principle, it is assumed that none of the elements will be divided. Then, the iterative process starts searching the k neighbors of every pressure element with a positive associated contact pressure (step B1). For this purpose, the algorithm proposed by Samet [19] is used, which is based in the quadtree data structure. This algorithm has been conveniently modified to account for several quadtree structures joined together. By using the quadtree data structure, the neighbors of a pressure element can be found although they do not share a vertex. As a result, this algorithm provides an array that contains the indexes of the k neighbors of a given pressure element.
The rate of change of the contact pressure u j;i between a pressure element D i and any of his neighbors D j is obtained using Eq. 8. If u j;i is lower than a user defined value u max , the next neighbor pressure element D jþ1 is evaluated. In contrast, if the rate of change of the contact pressure between both elements is greater than u max , these elements are candidates to be split. Three different situations arise: 1. Both pressure elements have the same level (L i ¼ L j ) and its value is lower than L max . Then, both elements are marked to be split, defining Both pressure elements have the same level (L i ¼ L j ) and its value is greater or equal to L max . Then, none of these elements are marked to be split, and the next neighbor is evaluated. 3. The pressure elements have different level (L i 6 ¼ L j ). If L i \L j the pressure element D i is marked to be split, defining K i ¼ TRUE. Otherwise, the pressure element D j is marked to be split, defining K j ¼ TRUE.
By using this criteria it is ensured that the level difference between two adjacent elements does not differ more than one level in the quadtree structure, avoiding unbalanced meshes.
The algorithm finishes when all the elements have been evaluated, returning an array which contains the indices of those elements where K i ¼ TRUE.

Final remarks
The topology of the resulting pressure element mesh, inside and outside the true contact area, depends on the configuration of the proposed approach, which is defined by a unique combination of the three input parameters: • The initial level of uniform mesh density, L ini • The maximum degree of mesh refinement, L max • The maximum allowed rate of change of contact pressure between adjacent elements, u max The possible configurations of the approach, and their effect on the resulting pressure element mesh, are categorized into three different settings, that are shown in Table 1. Strictly speaking, L ini must be greater than zero, although it is recommended that L ini ! 2, so a difference between neighbor pressure elements is guaranteed for the first iteration of the algorithm, even in axisymmetric contact problems.
Regarding the computational cost of the proposed approach, it has been observed that the most timeconsuming steps of the main algorithm are step A5 and step A6.
It can be demonstrated that the computer time invested in step A5 ðt A5 Þ is proportional to N f (Eq. 7), that is a quadratic function of the number of pressure elements involved in the solution of the contact problem. On the other hand, the computer time invested in step A6 ðt A6 Þ can be represented by a cubic function of the same number of elements. Taking this into account, it could be expected that t A6 [ t A5 when the number of pressure elements is relatively high. However, as it is explained by Kalker [11], the higher polynomial coefficients of t A5 compared to the ones of t A6 make t A5 much higher than t A6 for the vast majority of real cases.
In fact, in all the studied cases, t A5 represents between 80 and 95% of the computer time of the approach, in such a way that the former almost defines de latter. For this reason, and because t A5 is strictly related to N f , the parameter N f is considered a good descriptor of the computational cost of the algorithm.

Numerical examples and discussion
The performance of the proposed approach, in terms of accuracy and computational cost, is illustrated with two different cases of study. The case of study I corresponds to a punctual contact between a plane and a spheroidal indenter. The case of study II corresponds to a line contact between a plane and a cylindrical indenter. The dimensions of the indenters are shown in Fig. 6. A total contact load F T ¼ 40 kN is considered. The material of the plane is assumed to have a Young modulus of 70 GPa and a Poisson coefficient of 0.35. The material of both indenters (cases I and II) have a Young modulus of 210 GPa and a Poisson coefficient of 0.3.
In both cases, the root cell of the quadtree results in a 20 mm Â 20 mm square. The spheroidal indenter has been considered as an elastic half space. In contrast, finite dimensions have been considered for the longitudinal direction of the cylindrical indenter.
The accuracy and computational cost of the proposed approach are evaluated by solving the contact problems defined by cases of study I and II under several configurations of the approach, defined by its input parameters, that are shown in Table 2. For the aim of clarity, these configurations are categorized into several subgroups, where one of the input parameters is varied in a given range of values. The subgroups are categorized into three groups, where two input parameters can be varied. Each group allows for the study of the performance of the approach under one of the settings described in Table 1: • Group 1. A uniform mesh is used for the whole domain of the contact problem (Sect. 5.1). • Group 2. Adaptive mesh refinement is performed outside the true contact area (Sect. 5.2). • Group 3. Adaptive mesh refinement is performed both inside and outside the true contact area (Sect.

5.3).
In addition, the evolution of the computational cost with the size of the true contact area is also studied (Sect. 5.4). The computational cost of the approach to solve the contact problem under each configuration defined in Table 2 is measured using Eq. 7. On the other hand, the accuracy of the obtained solution is evaluated through the following parameters, that are measured as relative errors from the reference results shown in Table 3: • Relative error in maximum contact pressure, e R ðp max Þ • Relative error in the size of the true contact area, e R ðAÞ • Relative error in approximation between bodies, e R ðdÞ For case of study I, the reference results are determined using the analytical solution provided by the Hertz theory [10]. In constrast, Hertz theory is not longer applicable for case of study II, since its assumptions are not fulfilled because of the finite dimensions of the indenter in the longitudinal direction of the cylinder. For this reason, in this case the (I) (II)   In addition, in those cases where Hertz theory provides an analytical description of the contact pressure distribution (e.g. case of study I), an additional parameter is studied, that measures the error committed by the algorithm in the description of the contact pressure distribution. This parameter is defined as: where r 0 is the position vector of dA, and pðr 0 Þ is the theoretical contact pressure at position r 0 , determined from the Hertz theory.
5.1 Performance of the approach when a uniform mesh is used for the whole domain of the contact problem When a uniform pressure element mesh is used for the whole domain of the contact problem, a convergence of the results obtained by the proposed approach towards the reference results (Table 3) is expected as the pressure element mesh is refined (by increasing the value selected for L ini ). The performance of the approach under this situation is evaluated in this section. To do so, the contact problems defined by cases of study I and II are solved using the group 1 of configurations shown in Table 2. An example of the resulting contact area and pressure element mesh is shown in Fig. 7a.
The evolution of the parameters that define the accuracy of the approach as L ini is varied is shown in Fig. 8. As expected, it can be observed that as the pressure element mesh is refined, the results obtained by the proposed approach converge to the reference solution (i.e. the relative error tends to zero).
On the other hand, when a uniform mesh is used, the number of pressure elements in the discretization of the contact problem is n ¼ 4 L ini , and the computational cost of the algorithm is proportional to n 2 . For any value of L ini , the computational cost of the approach to solve the case of study II is always greater than the computational cost of the approach to solve the case of study I. This is because the factor of proportionality of the computational cost is 1 when no finite dimensions are taken into account, and 3 when two finite dimensions are taken into account, as happens in case of study II.

Performance of the approach when adaptive
mesh refinement is performed outside the true contact area In this section, the performance of the proposed approach when adaptive refinement is performed outside the true contact area is studied. To do so, the contact problems defined by cases of study I and II have been solved using the group 2 of configurations, shown in Table 2. Figure 7b shows an example of the resulting contact area and pressure element mesh when these configurations are used. The visual comparison between Fig. 7a and b allows to understand the reduction of pressure elements outside the true contact area under these configurations. The results obtained in these cases show that the accuracy in which the contact problem is solved is independent of the value selected for L ini . For any given value of L max , the same results as the ones shown in Fig. 8 (obtained with a uniform pressure element mesh for the whole domain of the contact problem) have been obtained, regardless of the value selected for L ini . This implies that the variations of the pressure element mesh outside the true contact area does not have any impact on the solution of the contact problem.
On the other hand, the evolution of the computational cost of the proposed approach to solve case of study I is shown in Fig. 9a. It can be observed that an important reduction of the computational cost is achieved as the difference between L max and L ini is increased, specially for those cases where the value selected for L max is large.
Although specifying a low value for L ini implies that more subdivisions of the pressure elements are    required to reach L max in all the elements of the true contact area, the number of influence coefficients that need to be calculated during the refinement is globally reduced, and so it is the computational cost of the proposed approach. However, the observed reduction of the computational cost is asymptotic, and once a certain value of L ini is reached, reducing L ini does not necessarily bring a significant reduction of the computational cost. Similar tendencies are observed for case of study II, although they are not shown for the aim of brevity.
In conclusion, it can be said that when u max ¼ 0, the computational cost of the proposed approach can be reduced by maximizing L max À L ini , while its accuracy remains unaffected.

Performance of the approach when adaptive
mesh refinement is performed both inside and outside the true contact area.
In this section, the performance of the proposed approach when adaptive refinement is also performed inside the true contact area is studied. To do so, the contact problems defined by cases of study I and II are solved using the group 3 of configurations shown in Table 2. An example of the resulting pressure element mesh is shown in Fig. 7c.
The obtained results show that the accuracy of the approach to predict the size of the true contact area does not depend on the value selected for u max , since the same values are obtained regardless of the value selected for this parameter. This is because the accuracy in which the border of the contact area is computed depends only on the value selected for L max , as stated in Sect. 4.1. In a similar way, it has been observed that the variation of the value selected for u max has a minor impact in the resulting approximation between bodies d, since the variation of the relative error e R ðdÞ with u max in any subgroup of configurations is lower than 0:01%.
On the other hand, Fig. 10a shows that the variation of u max produces changes in the accuracy in which the maximum contact pressure is calculated, since it depends on the resulting mesh distribution inside the contact area. However, as the value specified for u max gets lower, the result obtained for maximum contact pressure with the presented algorithm converges towards the reference solution.
In addition, Fig. 10b shows that the error committed in the description of the contact pressure distribution is slightly increased as it does the value selected for u max , specially in those cases where L max is large.
The consequences of the variations of u max over the contact pressure distribution are better observed in Fig. 11, where the calculated contact pressure distribution along the major and minor semiaxes of the contact ellipsis of case of study I are shown for two representative configurations of the approach, where only u max is varied. These contact pressure distributions correspond to the contact solutions shown in Fig  7. It can be observed that increasing the value selected for u max implies that a coarser mesh is used in those areas where the contact pressure gradient is small, without a significative loss of accuracy when describing the contact pressure distribution. However, in those cases where the maximum contact pressure is produced in an area where the pressure gradient is small, an increase of e R ðp max Þ can be expected.
Finally, Fig. 9b shows the evolution of the computational cost of the algorithm to solve the contact problem as u max is varied. It can be observed that the computational cost can be reduced by increasing the value selected for u max . Although this reduction is not as important as the one achieved by maximizing L max À L ini (Fig. 9a), it still can help to reduce the computational cost of the approach, specially in those cases where large values are selected for L max .
The results shown in this section show that the variation of u max has no effects in the solution of the contact problem when any configuration of subgroup G3.1 is used. This is because in this case, in all pressure elements within the contact area, L max is reached before the target value for u max is achieved.
In conclusion it can be said that specifying a value of u max [ 0 can help to reduce the computational cost of the approach. The accuracy in which the contact area is determined is unaffected by the variations of u max , but a loose of accuracy can be expected for the rest of the observed parameters.

Evolution of the computational cost
with the size of the true contact area When a uniform pressure element mesh is used to solve a contact problem with the proposed approach, its computational cost only depends of the value selected for L ini . In contrast, when adaptive refinement of the pressure element mesh is used, the computational cost is also dependent on the relation between the size of the true contact area (A) and the total size of the discretized domain ðA T Þ.
In this section, the evolution of the computational cost with the ratio A=A T is studied. For such a purpose, the case of study I has been solved varying the nominal value of the total contact force F T , and consequently, the size of the resulting true contact area is also varied, as shown in Fig. 12.
In first place, the evolution of the computational cost when adaptive refinement is only performed outside the true contact area is studied. Figure 13a shows the evolution of the computational cost when  Fig. 12 Axisymmetric representation of the resulting contact area for case of study I under several contact loads F T . These results are obtained using L ini ¼ L max ¼ 8 the contact problem is solved using the subgroup of configurations G2.4 described in Table 2. It can be observed that the benefit of using adaptive refinement outside the true contact area, in terms of computational cost, is greater as F T decreases (the ratio A=A T decreases).
In second place, the evolution of the computational cost when adaptive refinement is performed both inside and outside the true contact area is studied. Figure 13b shows the evolution of the computational cost of the approach when the contact problem is solved under the subgroup of configurations G3.4 described in Table 2. In these cases it can be seen that the benefit of using an adaptive refinement is greater as F T increases (the ratio A=A T increases).
These results lead to the conclusion that the benefit of performing adaptive refinement outside the true contact area is increased as its relative size is decreased. In addition, it has been observed that the benefit of performing adaptive refinement inside the true contact area is increased as its relative size is increased.

Conclusions
In this work, a new approach has been developed to solve non-hertzian contact problems through the influence coefficient method, that improves the classical approach by using an adaptive mesh refinement that is based in a quadtree decomposition of the domain. Starting from an initial level of uniform mesh density (defined by the parameter L ini ), a mesh refinement is performed based on two different criteria: (1) the maximum allowed rate of change of contact pressure between adjacent elements, defined by the parameter u max , and (2) the maximum degree of mesh refinement, defined by the parameter L max .
The configuration of the approach is defined by a unique combination of values for L ini , L max and u max . Depending on its configuration, the approach performs a mesh refinement, that can be performed only outside the true contact area, or both inside and outside the true contact area.
The proposed approach has been tested through several study cases and configurations. The obtained results, in terms of maximum contact pressure, true contact area and approximation between bodies, enable us to draw the following conclusions: • When L ini ¼ L max , a uniform mesh is used to solve the contact problem, regardless of the value selected for u max . Under this configuration, it can be observed that the obtained results converge towards the reference solution as L ini is increased. However, the computational cost increases exponentially. • When L ini \L max and u max ¼ 0, adaptive mesh refinement is performed outside the true contact area. Under these circumstances, it can be observed that the computational cost of the approach is reduced by maximizing the ratio L max À L min , while the accuracy of the solution remains unaffected. • In last place, when L ini \L max and u max [ 0, adaptive mesh refinement is performed both inside and outside the true contact area. Under these circumstances, it can be observed that a further reduction of the computational cost can be achieved. However, although the accuracy in the predicition of the size of the contact area and  Fig. 13 Influence of the size of the true contact area on the computational cost of the proposed approach approximation between bodies is not affected by u max , a loss of accuracy can be expected in the prediction of the contact pressure distribution as u max is increased.
In terms of computational cost, the benefit of using adaptive mesh refinement outside the true contact area is greater as the size of the true contact area decreases.
On the other hand, the benefit of using adaptive refinement inside the true contact area increases as the size of the true contact area increases. In general, it can be said that when using the proposed approach to solve frictionless elastic nonhertzian contacts, maximizing L max À L ini and specifying a low value for u max can yield an important reduction of the computational cost without a significant loss of accuracy.
Further works on this topic could go directed towards adapting the proposed mesh refinement strategy to solve contact problems in which the friction effect and the tangential behavior need to be considered. Apart from that, the investigation of other refinement criteria or the exploration of other adaptive mesh refinement strategies (such as the anisotropic mesh adaptation) could be an interesting extension of this work.