Open access peer-reviewed chapter

# Adaptive Mesh Refinement Strategy for the Semi-Analytical Solution of Frictionless Elastic Contact Problems

Written By

Victor Roda-Casanova and Francisco Sanchez-Marin

Submitted: 21 April 2017 Reviewed: 13 November 2017 Published: 20 December 2017

DOI: 10.5772/intechopen.72422

From the Edited Volume

## Contact and Fracture Mechanics

Edited by Pranav H. Darji and Veera P. Darji

Chapter metrics overview

View Full Metrics

## Abstract

Semi-analytical methods are commonly used to solve contact problems. These methods require the discretization of the domain into a mesh of pressure elements. In general, it can be said that their accuracy increases as the pressure element mesh is refined. However, the refinement of the pressure element mesh also implies an increase in their computational cost. So, in the great majority of the cases, a commitment between accuracy and computational cost must be achieved. In this work, a new approach is presented, whose main purpose is to improve the efficiency of the semi-analytical methods that are used to solve contact problems. To do so, an adaptive refinement of the pressure element mesh is implemented. This strategy allows for a reduction of the computational cost of the method, while its accuracy remains unaffected.

### Keywords

• contact analysis
• semi-analytical methods

## 1. Introduction

The contact stress analysis plays an important role during the design process of several mechanical elements like bearings, gears, etc. In order to accomplish a contact analysis, the so-called contact problem must be solved to obtain the following relevant information:

1. The contact area, which involves the determination of the size, shape, and location of the true contact area in each one of the contacting bodies.

2. The contact stresses, which involve the determination of the contact pressure distribution on the surface of the bodies and the stress distribution underneath the surfaces.

3. The deformation of the bodies produced by the contact pressure.

Different approaches have been used to solve contact problems, which can be classified into three groups: numerical, analytical, and semi-analytical methods. Compared to the numerical methods, it can be said that the analytical methods are more efficient in terms of computational cost, but they have severe applicability limitations imposed by the hypotheses of the underlying theory. On the other hand, the numerical methods can overcome these limitations, but at a much higher computational cost.

The semi-analytical methods (SAMs) can be considered as an intermediate approach: they are potentially faster than the numerical methods, while they allow overcoming some of the limitations of the analytical methods. SAMs are usually based on the discretization of the potential contact area into a mesh of n pressure elements, with a uniform pressure distribution assumed to be acting over each one of them. Influence coefficients are used to relate the pressure applied over each pressure element with the displacements that this pressure produces at the centroid of the other elements of the mesh. Using these influence coefficients, the solution to the contact problem can be numerically found in terms of the contact pressure distribution that satisfies the contact conditions.

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 [1] stated that the computational cost of these semi-analytical methods can be defined by the number of influence coefficients that need to be calculated to solve the contact problem (that, in general, is proportional to n2). He also argued that the accuracy of the solution to the contact problem, in terms of contact area and contact pressure distribution, depends on the refinement of the pressure element mesh, especially in those regions 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 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 those cases, it is common to use a uniform pressure element mesh for the whole domain, being more or less dense depending on the desired accuracy and on the capabilities of the computer used to solve the contact problem. In consequence, there could be many pressure elements in the discretization out of the true contact area, what causes a loss 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 numerical methods based in the discretization of the domain, especially in FEM procedures [2]. However, no previous use of adaptive refinement has been found in the literature for the solution of contact problems using semi-analytical methods.

In this work, an approach to solve frictionless elastic contact problems is presented, whose main purpose is to improve the efficiency of the semi-analytical methods that are used to solve contact problems. To do so, an adaptive refinement of the pressure element mesh is implemented, which is based on the discrete rate of change of any magnitude that is related with the solution of the contact problem. This strategy allows for a reduction of the computational cost of the method, while its accuracy remains unaffected. The theoretical background and the computational implementation of the method are described, and its performance is illustrated with numerical examples.

## 2. Theoretical background

This section describes the theoretical background under which the proposed approach to solve frictionless elastic contact problems is developed. The concept of pressure element is described, as well as those considerations required to solve contact problems between bodies of finite dimensions. Finally, the quadtree decomposition of the domain is introduced, which is a useful strategy to perform adaptive mesh refinement.

### 2.1. Pressure elements and surface normal deflection in an elastic half-space

Consider a body that, because of its main features, can be approached to an elastic half-space, as the one shown in Figure 1a. A Cartesian coordinate system is defined over the surface of this body, which X and Y axes define a plane that is coincident with its surface, and the Z axis points inward him. A normal pressure distribution p is applied over the surface of the body, acting over an area that is denoted by S.

Now consider a generic point C within the area S, whose position is defined by the vector rxyz, being z=0. Consider another point H in the surface of the body, whose position is defined by the vector rxyz, being z=0. The normal elastic deflection produced at a point H due to a normal pressure distribution applied over the area S can be determined by the superposition of the Boussinesq relation [3]:

ωr=1ν2πGSprrrdSE1

where ν is the Poisson coefficient and G is the shear modulus of the material of the considered body.

Obtaining a generic closed-form solution for Eq. (1) is not possible, since it depends on the shape of the area S and on the considered pressure distribution. However, several closed-form solutions can be found in the literature for certain pressure distributions applied over areas with a specific shape (such as triangles, rectangles, hexagons, etc.).

Let’s focus on the closed-form solution for Eq. (1) that Love [4] obtained for uniform pressure distributions acting over areas with rectangular shape, as the one shown in Figure 1b. From now on, this combination of shape and pressure distribution will be called pressure element, and will be denoted by Δj. The area of the pressure element shown in Figure 1b is Aj=2a×2b, and the uniform pressure distribution that acts over this area is pr=pj. Under these conditions, the closed-form solution for Eq. (1) is

ωr=fjrpjE2

where fjr is the influence coefficient of pressure element j over the point H, which can be analytically determined as

fjr=1ν2πGClnA+A2+C2B+B2+D2+AlnC+A2+C2D+A2+D2+DlnB+B2+D2A+A2+D2+BlnD+B2+D2C+B2+C2E3

where coefficients A, B, C, and D are calculated as

A=dy+bC=dx+aB=dybD=dxa

These pressure elements can be useful to determine the normal displacement produced at the surface of a body due to a non-uniform pressure distribution applied over a complex area. To illustrate this methodology, consider a complex area S, as the one shown in Figure 2a, over which an arbitrary pressure distribution is acting. To determine the displacement field produced by this pressure distribution, the area S is discretized into a mesh of n rectangular pressure elements Δj, as shown in Figure 2b. Then, the arbitrary pressure distribution is approached by assigning a uniform pressure value pj to each pressure element, as shown in Figure 1c.

Finally, the displacement at any point of the surface of the body can be determined by superposition of the displacements produced at this point by uniform pressures acting over each pressure element of the mesh as

ωr=j=1npjfjrE4

This methodology can be applied with different types of pressure elements, having different shapes and pressure distributions acting over them. However, using rectangular pressure elements has some advantages, which are discussed in [5].

### 2.2. Semi-analytical method to solve frictionless elastic contact problems

The solutions exposed in the previous section can be used to obtain the pressure distribution that is produced when two bodies are pressed together in the absence of friction. For such a purpose, it is necessary that the two bodies can be approached to elastic half-spaces in the vicinity of the area in which the contact between them is produced.

Consider two bodies 1 and 2 in its undeformed contact position, contacting at the initial point of contact OL (Figure 3a). At this point, a common tangent plane Π is defined, which 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 can be referred to it in the linear small strain theory of elasticity.

A Cartesian coordinate system is defined with origin at point OL, being the local axis ZL normal to the plane Π and pointing inward the body 2. Consider a generic point Q in the plane Π, whose position is defined by the vector rxLyLzL, being zL=0. The gap between the two bodies, measured along ZL axis, is denoted by the function B r, which in the first instance is assumed to be smooth and continuous.

The two bodies are pressed together in the absence of friction by the effect of the force FT (Figure 3b), causing a normal approach between them that is denoted by δ. Since penetration is physically inadmissible, a contact pressure distribution pr is generated in the true contact area S that deforms the contacting bodies. In this way, elastic normal deflections are produced in the surfaces of the bodies 1 and 2, which are denoted by ω1r and ω2r, respectively. The total normal deflection is denoted by ωr and can be calculated as

ωr=ω1r+ω2rE5

The essence of the resulting contact problem is to determine the pressure distribution that fulfills the contact conditions, both inside and outside the contact area, whose geometry and size are unknown

pr>0andBr+ωrδ=0insideSpr=0andBr+ωrδ>0outsideSE6

Kalker [5] demonstrated that the solution to this contact problem can be found minimizing the total complementary energy V under the condition that the contact pressure is equal or greater than zero in all the domain of the problem. The total complementary energy is defined as the sum of the internal complementary energy of the stressed bodies and the external complementary energy as

V=12SprωrdS+SprBrδdSE7

To enable the numerical solution, the potential contact area is discretized into a set of n pressure elements Δj (described in Section 2.1), with a uniform pressure distribution assumed to be acting over each one of them, as shown in Figure 4. The position of the centroid of each pressure element Δj is denoted by vector rj.

Under a discretized domain, the total complementary energy may be expressed as

V=12i=1npiAiωrdAi+i=1npiAiBrδdAiE8

Taking into account Eq. (4), Eq. (5) may be rewritten as

ωr=j=1npjfj1r+j=1npjfj2r=j=1npjtjrE9

where tjr is defined as the cumulated influence coefficient of pressure element Δj over point Q

tjr=fj1r+fj2rE10

Considering Eq. (9), Eq. (8) can be rewritten as:

V=12i=1npij=1npjAitjrdAi+i=1npiAiBrδdAiE11

To reduce the computational cost of the calculations, two assumptions are made:

1. The distance between the surfaces of the bodies Br is assumed to be constant in all the pressure element Δi, and equal to the distance between the surface of the two bodies at its centroid, in such a way that Br=Bri=Bi.

2. The cumulated influence coefficient tjr is assumed to be constant over all the pressure element Δi, and equal to the value in its centroid, in such a way that tjr=tjri=tj,i. The coefficient tj,i can be defined as the cumulated influence coefficient of element Δj over the centroid of element Δi, and it may be expressed as:

tj,i=fj,i1+fj,i2

where fj,i1 and fj,i2 are the influence coefficients of element Δj over the centroid of element Δi, which can be determined for each contacting body using Eq. (3).

Under these assumptions, the total complementary energy can be expressed as

V=12i=1nj=1npipjtj,iAi+i=1npiBiδAiE12

The solution to the contact problem, in terms of contact pressure distribution, can be found by minimizing Eq. (12) under the following restrictions:

Vpi=j=1npjtj,iAi+BiδAi=0ifpi>0Vpi=j=1npjtj,iAi+BiδAi0ifpi=0E13

The true contact area is then defined, within the precision of the mesh size, by the boundary between the elements with zero and non-zero pressures. The total contact load FT can be calculated as

FT=i=1npiAiE14

### 2.3. Contact between bodies of finite dimensions

The method described in the previous section is based on the utilization of influence coefficients that relate the contact pressures with the surface displacements. These influence coefficients can be determined in several ways, but in the great majority of the cases, they are calculated using the superposition of the Boussinesq relation (described in Section 2.1).

Because of the hypotheses under which this relation is established, influence coefficients determined using the Boussinesq relation should only be used to solve contact problems between contact bodies that can be approached to half-spaces. Otherwise, the application of the described method can lead to erroneous solutions of the contact problem.

However, the influence coefficients determined using the Boussinesq relation can be corrected, so they can be used to solve contact problems between bodies that a priori cannot be approached to elastic half-spaces. Among the correction methods that can be found in the literature, a correction method to consider contact bodies of finite dimensions is described in this section.

In this correction method, the finiteness of the contacting bodies is characterized by stress-free surfaces that are perpendicular to the length direction of the bodies, as illustrated in Figure 5. To leave those areas of the half-space that coincide with these surfaces free of normal and shear stresses, the correction method proposed by de Mul [6] is used, which consist in modifying the calculation of the influence coefficients fj,i for each body with finite dimensions in the following way:

fj,i=fjo,i+1.2911ν0.080.5νΨm=1nfjm,iE15

where fj0,i is the influence coefficient of the original pressure element Δjo, Ψ is a correction factor proposed by Guilbault [7], and fjm,i are the influence coefficients of the pressure elements Δjm, that are mirrored instances of the element Δjo respect to the planes that coincide with the n free-stress surfaces of the body.

This method involves the calculation of additional influence coefficients of the mirrored pressure elements, and hence the computational cost of the method is multiplied by n+1, being n the number of finite dimensions taken into account in the problem.

### 2.4. Quadtree decomposition of the domain

According to Samet [8], the basic concept of the quadtree is to enclose the domain of the problem Γ into a containing cell, usually a square, which is denoted as the root of the quadtree, as shown in Figure 6a. This cell is then subdivided into four sons of the same size (Figure 6b), 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 (Figure 6c and d).

The information related to the quadtree decomposition of the domain is stored in a hierarchical tree structure, as shown in Figure 6e. For every cell, references to its ancestor and sons are stored. This kind of structure eases the performance of several operations, such as the neighbor finding in a defined direction, which will play 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 Lj, and represents the number of divisions performed from the root of the quadtree. According to this definition, Lj 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 Lj 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 gray in Figure 6e), while subdivided cells are referred to as non-leaf cells.

## 3. Computational approach to solve frictionless elastic contact problems using adaptive mesh refinement

In this section, the major topics of the computational implementation of the semi-analytical method to solve frictionless elastic contact problems described in Section 2.2 are discussed. As it has been said before, this method is based on the discretization of the potential contact area Γ into a mesh of pressure elements. The potential contact area is defined as the Boolean intersection of the projection of the bodies on the plane Π, as shown in Figure 7.

In the classical approach, the potential contact area is discretized using a uniform mesh of pressure elements. In contrast, to improve the efficiency of the method, adaptive mesh refinement is implemented in this approach. To do so, a quadtree decomposition (described in Section 2.4) of the potential contact area is performed, where all the leaf cells of the quadtree are considered pressure elements. The use of a quadtree offers two interesting features to this implementation. In first place, the recursive division of the cells provides a robust local mesh refinement strategy. In second place, transverse operations such as neighbor finding algorithms are computationally efficient and easy to implement.

The main algorithm of the approach to solve contact problems using adaptive mesh refinement is shown in Figure 8. The following inputs are required by the algorithm:

1. The geometry and position of the contact surfaces in undeformed contact position.

2. The initial point of contact OL and a vector defining the contact normal.

3. The magnitude of the contact force FT.

4. The initial level of uniform mesh density Luni, which 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 Π, where a local Cartesian coordinate system is defined, being the ZL axis normal to the plane Π (step A1). The boundaries of the contacting bodies are normally projected onto the plane Π to determine the potential contact area Γ (step A2).

The potential contact area is enclosed by the root cell of a quadtree, which is recursively subdivided until the desired initial level of uniform mesh density Luni is reached for all the cells of the quadtree (step A3). All leaf cells of the quadtree are considered pressure elements Δi for the initial iteration of the algorithm, which is performed using a uniform pressure element mesh. All the elements are marked with the flag Λi=TRUE, indicating that their properties (associated area Ai, normal gap Bi, contact pressure pi, and cumulated influence coefficients tj,i) are not computed yet.

To maximize the efficiency of the proposed approach, it is important to minimize the number of pressure elements located outside of the potential contact area. This can be achieved by ensuring that the potential contact 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 [9].

Then, an iterative process starts whose first step is the determination of the normal gap Bi for all the pressure elements (step A4) where Λi=TRUE. The cumulated influence coefficients tj,i of those elements are also determined (step A5), using Eq. (3). Finally, the contact problem is solved using Eq. (13) (step A6), in the way indicated in [5], obtaining a contact pressure value pi for each element Δi of the discretization.

The flag Λi is defined as FALSE for all the elements present in the discretization, indicating that the properties of these elements have already been computed.

At this point, the adaptive mesh refinement is performed. For such a purpose, the algorithm that determines the elements that must be split (that is described in section 3.1) is called (step A7), which returns an array with the indexes of these elements. Then, the selected 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 Λ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 implementation 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 method. The number of influence coefficients calculated in the proposed approach Nf can be determined afterwards using the following equation:

Nf=i=1t2ninnewinnewi2E16

where t is the number of iterations performed by the algorithm, ni is the number of elements in iteration i, and nnewi is the number of new elements in iteration i.

### 3.1. Algorithm to determine elements to split

An adaptive mesh refinement may be based upon several criteria. In this work, the rate of change of a given physical magnitude (denoted by λ) related with the solution of the contact problem is used to perform adaptive refinement of the pressure element mesh. Since the proposed approach works under a discretized domain, each pressure element Δi will have an associated value λi of the observed physical magnitude, and in consequence, a discrete rate of change φj,i of λ can be established between an element Δi and any of its neighbors Δj as

φj,i=λjλimaxλjλiE17

If the discrete rate of change φj,i between a pressure element Δi and any of its neighbors Δj is higher than an arbitrarily defined value φmax (representing the maximum allowed rate of change of the observed physical magnitude), then the pressure element Δi is marked as a candidate to be split.

However, in some situations, the rate of change of λ is so high that the condition φj,i<φmax cannot be reached, and the refinement strategy based on the rate of change of λ would refine the pressure element mesh endlessly. In order to limit the number of iterations performed by the algorithm, an additional stopping criterion based on the minimum size allowed for a pressure element needs to be included. As mentioned before, the level Lj 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 Lmax, referred to as the maximum degree of mesh refinement.

It is important to point out that φmax is a target value and may not be always reached. If the maximum degree of mesh refinement Lmax is reached for all pressure elements before the target value for φmax is achieved, the mesh refinement will finish.

The main routine of the algorithm to determine elements to split is shown in Figure 8b. The following input information is required by the algorithm:

1. The properties associated with the pressure elements (area Ai, normal gap Bi, contact pressure pi, and cumulated influence coefficients tj,i).

3. The maximum degree of mesh refinement Lmax.

4. The maximum allowed rate of change of the observed physical magnitude φmax.

The algorithm starts defining the flag Ki as FALSE for all the elements present in the current discretization. The flag Ki 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 Δi in the domain (step B1). For this purpose, the algorithm proposed by Samet [8] is used, which is based in the quadtree data structure. As a result, this algorithm provides an array that contains the indexes of the k neighbors of a given pressure element.

For each pair of neighboring pressure elements Δi and Δj, the algorithm determines the associated physical magnitudes λi and λj (step B2). These magnitudes can be already given by the main algorithm (as in the case of the contact pressures), or they can be specifically determined from these values by performing additional calculations.

Then, the discrete rate of change φj,i of the observed magnitude between pressure element Δi and his neighbor Δj is obtained using Eq. (17) (step B3). If φj,i is lower than the user-defined value φmax, the next neighbor pressure element Δj+1 is evaluated. In contrast, if the discrete rate of change of the observed magnitude between both elements is greater than φmax, then element Δi is considered as a candidate to be split. Two additional conditions must be fulfilled so Δi can be marked to be split:

1. On one hand, Li must be lower than Lmax, to avoid that the algorithm refines the mesh indefinitely in those cases where φmax cannot be reached.

2. On the other hand, Li must be less than or equal to Lj, to ensure that the level difference between two adjacent elements does not differ more than one level in the quadtree structure, avoiding unbalanced meshes.

If these conditions are met, the pressure element Δi is marked to be split by defining Ki=TRUE. The algorithm finishes when all the pressure elements have been evaluated, returning an array that contains the indices of those elements where Ki=TRUE.

### 3.2. Final remarks

In contact problems, there are different physical magnitudes that can be observed to perform the adaptive mesh refinement, having each one its advantages and disadvantages. In this work, the observed magnitude to perform the adaptive mesh refinement is the contact pressure (λi=pi), and the mesh refinement is performed based on the discrete gradient of the contact pressures. The main advantages of choosing the gradient of the contact pressure as refinement criterion instead of any other derived magnitude are:

1. On one hand, in this approach, the solution of the contact problem is found in terms of the contact pressure distribution. From the calculated contact pressure distribution, derived results are obtained. Since the accuracy of the derived results is dependent from the accuracy in which contact pressure distribution is calculated, it is important to obtain an accurate description of the contact pressure distribution.

2. On the other hand, using the contact pressure distribution, instead of the derived results, as refinement criteria helps reducing the computational cost of the proposed approach, because obtaining derived results implies additional calculations.

However, it must be taken into account that the contact pressure distribution function is not differentiable in the border of the contact area. In consequence, according to Eq. (17), the discrete rate of change of the contact pressure between an element Δi that is within the contact area (pi>0) and of an adjacent element Δj that is outside of the contact area (pj=0) is always φj,i=1. Therefore, if a value lower than 1 is specified for φmax, the refinement strategy will refine the mesh at the boundary of the contact area until the maximum degree of mesh refinement will be reached at the border of the true contact area.

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:

1. The initial level of uniform mesh density, Luni.

2. The maximum degree of mesh refinement, Lmax.

3. The maximum allowed rate of change of the physical magnitude, φmax.

The possible configurations of the approach, and their effect on the resulting pressure element mesh, are categorized intro three different settings:

1. Setting 1 (LuniLmax, φmax not relevant): using this setting, the contact problem is solved using a uniform mesh, whose mesh density is defined by Luni.

2. Setting 2 (Luni<Lmax,φmax=0): using this setting, the contact problem is solved using adaptive mesh refinement outside the true contact area. Inside the true contact area, a uniform mesh is used, whose mesh density is defined by Lmax.

3. Setting 3 (Luni<Lmax,φmax>0): using this setting, the contact problem is solved using adaptive mesh refinement both inside and outside the true contact area.

From the nine steps of the main algorithm, step A5 is the most time consuming. For this reason, the computational cost of the approach can be defined by the number of influence coefficients that are calculated to solve the contact problem, which can be determined using Eq. (16).

## 4. Numerical examples

The performance of the proposed approach is illustrated in this section, considering its accuracy and computational cost. For such a purpose, two cases of study are considered:

• Case of study I (CoSI) corresponds to a punctual contact between a plane and a spherical indenter, whose dimensions are shown in Figure 9a. Punctual contacts are common in mechanical components such as ball bearings, gears and rail-wheel systems.

• Case of study II (CoSII) corresponds to a line contact between a plane and a cylindrical indenter, whose dimensions are shown in Figure 9b. Line contacts are common in mechanical, such as roller bearings or standard spur and helical gears.

The material of both indenters (CoSI and CoSII) and the plane is assumed to have a Young modulus of 70GPa and a Poisson coefficient of 0.35. A total contact load FT=60kN is considered.

In both cases, the root cell of the quadtree results in a 20×20mm square. The spherical indenter has been considered as an elastic half-space. In contrast, two finite dimensions have been considered for the longitudinal direction of the cylindrical indenter, using the correction method described in Section 2.3.

The cases of study I and II are solved under several configurations of the proposed approach, selected from the three settings described in Section 3.2, and the performance of each configuration is discussed in Sections 4.1 (for configurations within setting 1), 4.2 (for configurations within setting 2), and 4.3 (for configurations within setting 3).

For each configuration, the computational cost of the approach to solve the contact problem is evaluated using Eq. (16). The accuracy of the approach is evaluated by comparing the obtained contact pressure distributions with reference solutions. For case of study I, the reference solution is determined using the analytical solution provided by the Hertz contact theory [10]. In contrast, since Hertz theory is no longer applicable for case of study II, reference results are obtained for this case using a validated finite element model.

### 4.1. Performance of the approach when a uniform mesh is used for the whole domain of the contact problem

The performance of the approach when a uniform pressure element mesh is used for the whole potential contact area is illustrated in this section. To do so, the contact problems defined by cases of study I and II are solved under several configurations of the approach, in which Luni has been varied, keeping Luni=Lmax and φmax=0 (setting 1 in Section 3.1). Figure 10ac show examples of the resulting contact area and pressure element mesh that have been obtained for case of study I under this setting of the approach. The computational cost of the proposed approach to solve the case of study I is also shown for each configuration.

The contact pressure distributions along the principal axes of the contact area of the solutions shown in Figure 10ac are shown in Figure 11a. As expected, it can be observed that as the pressure element mesh is refined (by increasing the value selected for Luni), the results obtained by the proposed approach converge toward the reference solution.

Using this configuration of the approach, a mesh containing 4Luni pressure elements is used, regardless of the nature of the contact problem to be solved. Under these circumstances, the computational cost is proportional to 42Luni, and the factor of proportionality is the number of finite dimensions taken into account in the contact problem (as explained in Section 2.3). In consequence, for any value of Luni, the computational cost of the algorithm to solve case of study II will always be greater than the computational cost to solve case of study I.

### 4.2. 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 illustrated. To do so, the contact problems defined by cases of study I and II are solved under several configurations of the approach, in which Luni and Lmax have been varied, keeping Luni<Lmax and φmax=0 (setting 2 in Section 3.1). Figures 10d and 12a show examples of the resulting contact area and pressure element mesh that have been obtained for cases of study I and II under this setting of the approach.

The results obtained in these cases show that the accuracy in which the contact problem is solved is independent of the value selected for Luni. For any given value of Lmax, the same contact pressure distributions as the ones obtained with a uniform pressure element mesh for the whole domain have been obtained (shown in Figure 11a), regardless of the value selected for Luni. This implies that the variation 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, comparing the computational cost of the solutions shown in Figure 10c (uniform mesh) and 10d (adaptive refinement outside the true contact area), it can be observed that an important reduction of the computational cost is achieved by increasing the difference between Lmax and Luni. Similar tendencies are observed for case of study II, where the reductions of computational cost are even more remarkable due to the presence of finite dimensions.

### 4.3. 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 illustrated. To do so, the contact problems defined by cases of study I and II are solved under several configurations of the approach, in which φmax has been varied, keeping Luni<Lmax (setting 3 in Section 3.1). Figures 10e, f and 12b show examples of the resulting contact area and pressure element mesh that have been obtained for cases of study I and II under this setting of the approach.

The contact pressure distributions along the principal axes of the contact area of the solutions shown in Figure 10df are shown in Figure 11b. The contact pressure distributions along the principal axes of the contact area of the solutions shown in Figure 12a, b are shown in Figure 13. It both cases, it can be observed that increasing the value selected for φmax implies that a coarser mesh is used in those regions of the true contact area where the contact pressure gradient is small, without a significant loss of accuracy when describing the contact pressure distribution.

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 φmax, since the same values are obtained regardless of the value selected for this parameter. This is because when φmax<1, the accuracy in which the border of the contact area is computed depends only on the value selected for Lmax, as stated in Section 3.2.

Finally, comparing the computational cost of the solutions shown in Figure 10d and e (and Figure 12a and b), it can be observed that a further reduction of the computational cost can be achieved by specifying values of φmax>0. Although this reduction is not as important as the one achieved by maximizing LmaxLuni (discussed in Section 4.2), it still can help to reduce the computational cost of the approach.

## 5. Conclusions

A new semi-analytical approach has been developed to solve frictionless elastic contact problems using adaptive mesh refinement. Starting from a coarse initial uniform mesh (whose density is defined by the parameter Luni), a mesh refinement is performed based on two different criteria: (i) the maximum allowed rate of change of a physical magnitude (the contact pressure), defined by the parameter φmax and (ii) the maximum degree of mesh refinement, defined by the parameter Lmax.

The configuration of the approach is defined by a unique combination of values for Luni, Lmax, and φmax. The performance of the proposed approach has been illustrated with several cases of study solved under different configurations of the approach, and the obtained results enable us to draw the following conclusions:

1. When Luni=Lmax, a uniform mesh is used to solve the contact problem, regardless of the value selected for φmax. Under this configuration, it can be observed that the obtained results converge toward the reference solution as Luni is increased. However, an exponential growth of the computational cost is produced as the pressure element mesh is refined.

2. When Luni<Lmax and φ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 LmaxLini, while the accuracy of the solution remains unaffected.

3. In last place, when Luni<Lmax and φ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, a loss of accuracy can be expected in the prediction of the contact pressure distribution as φmax is increased.

A further discussion on this topic can be found in Ref. [11].

## Acknowledgments

Parts of this chapter are reproduced from authors’ previous publication (Ref. [11]).

## References

1. 1. Kalker JJ, Van Randen Y. A minimum principle for frictionless elastic contact with application to non-Hertzian half-space contact problems. Journal of Engineering Mathematics. 1972;6(2):193-206
2. 2. Wriggers P. Computational Contact Mechanics. 2nd ed. Berlin/Heidelberg: Springer-Verlag; 2006. p. 518
3. 3. Boussinesq J. Application des potentiels à l’étude de l’équilibre et du mouvement des solides élastiques. Paris: Gauthier-Villars; 1885
4. 4. Love AEH. A Treatise on the Mathematical Theory of Elasticity. Cambridge: Cambridge University Press; 1906
5. 5. Kalker JJ. Three-Dimensional Elastic Bodies in Rolling Contact. 1st ed. Netherlands: Springer; 1990. p. 314
6. 6. de Mul JM, Kalker JJ, Fredriksson B. The contact between arbitrarily curved bodies of finite dimensions. Journal of Tribology. 1986;108(1):140-148
7. 7. Guilbault R. A fast correction for elastic quarter-space applied to 3D modeling of edge contact problems. Journal of Tribology. 2011;133(3):031402
8. 8. Samet H. Neighbor finding techniques for images represented by quadtrees. Computer Graphics and Image Processing. 1982;18:37-57
9. 9. Chaudhuri D, Samal A. A simple method for fitting of bounding rectangle to closed regions. Pattern Recognition. 2007;40(7):1981-1989
10. 10. Johnson KL. Contact Mechanics. Cambridge/New York: Cambridge University Press; 1985
11. 11. Roda-Casanova V, Sanchez-Marin F. Meccanica. 2017. https://doi.org/10.1007/s11012-017-0806-y

Written By

Victor Roda-Casanova and Francisco Sanchez-Marin

Submitted: 21 April 2017 Reviewed: 13 November 2017 Published: 20 December 2017