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

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.

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 $n^2$). 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 $\left(p\right)$ is applied over the surface of the body, acting over an area that is denoted by $S$.

Figure 1.

Pressure distributions applied over an elastic half-space.

Now consider a generic point $C$ within the area $S$, whose position is defined by the vector ${\mathit{r}}^{\prime }\left(x,y,z\right)$, being $z=0$. Consider another point $H$ in the surface of the body, whose position is defined by the vector $\mathit{r}\left(x,y,z\right)$, 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]:

$$\omega \left(\mathit{r}\right)=\frac{1-\nu }{2\mathit{\pi G}}{\int }_S\frac{p\left({\mathit{r}}^{\prime }\right)}{\left|\mathit{r}-{\mathit{r}}^{\prime }\right|}\mathit{dS}$$E1

where $\nu$ 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 ${\mathrm{\Delta }}_j$. The area of the pressure element shown in Figure 1b is $A_j=2a\times 2b$, and the uniform pressure distribution that acts over this area is $p\left({\mathit{r}}^{\prime }\right)=p_j$. Under these conditions, the closed-form solution for Eq. (1) is

$$\omega \left(\mathit{r}\right)=f_j\left(\mathit{r}\right)\bullet p_j$$E2

where $f_j\left(\mathit{r}\right)$ is the influence coefficient of pressure element ${∆}_j$ over the point $H$, which can be analytically determined as

$$\begin{array}{ll}\hfill & f_j\left(\mathit{r}\right)=\frac{1-\nu }{2\mathit{\pi G}}\left\{C\bullet \mathit{ln}\left[\frac{A+\sqrt{A^2+C^2}}{B+\sqrt{B^2+D^2}}\right]+A\bullet \mathit{ln}\left[\frac{C+\sqrt{A^2+C^2}}{D+\sqrt{A^2+D^2}}\right]+D\bullet \mathit{ln}\left[\frac{B+\sqrt{B^2+D^2}}{A+\sqrt{A^2+D^2}}\right]+B\bullet \mathit{ln}\left[\frac{D+\sqrt{B^2+D^2}}{C+\sqrt{B^2+C^2}}\right]\right\}\end{array}$$E3

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

$$\begin{array}{cc}A=d_y+b& C=d_x+a\\ B=d_y-b& D=d_x-a\end{array}$$

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 ${\mathrm{\Delta }}_j$, as shown in Figure 2b. Then, the arbitrary pressure distribution is approached by assigning a uniform pressure value $p_j$ to each pressure element, as shown in Figure 1c.

Figure 2.

Normal deflection produced by a complex pressure distribution.

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

$$\omega \left(\mathit{r}\right)=\sum _{j=1}^n\left[p_j\bullet f_j\left(\mathit{r}\right)\right]$$E4

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 $O_L$ (Figure 3a). At this point, a common tangent plane $\Pi$ 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.

Figure 3.

Contact between two bodies: (a) undeformed position and (b) deformed position.

A Cartesian coordinate system is defined with origin at point $O_L$, being the local axis $Z_L$ normal to the plane $\Pi$ and pointing inward the body $2$. Consider a generic point $Q$ in the plane $\Pi$, whose position is defined by the vector $\mathit{r}\left(x_L,y_L,z_L\right)$, being $z_L=0$. The gap between the two bodies, measured along $Z_L$ axis, is denoted by the function B $\left(\mathit{r}\right)$, 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 $F_T$ (Figure 3b), causing a normal approach between them that is denoted by $\delta$. Since penetration is physically inadmissible, a contact pressure distribution $p\left(\mathit{r}\right)$ 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 ${\omega }^{\left(1\right)}\left(\mathit{r}\right)$ and ${\omega }^{\left(2\right)}\left(\mathit{r}\right)$, respectively. The total normal deflection is denoted by $\omega \left(\mathit{r}\right)$ and can be calculated as

$$\omega \left(\mathit{r}\right)={\omega }^{\left(1\right)}\left(\mathit{r}\right)+{\omega }^{\left(2\right)}\left(\mathit{r}\right)$$E5

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

$$\begin{array}{cc}p\left(\mathit{r}\right)>0\text{and}B\left(\mathit{r}\right)+\omega \left(\mathit{r}\right)-\delta =0& \begin{array}{cc}\text{inside}& S\end{array}\\ \begin{array}{ccc}p\left(\mathit{r}\right)=0& \text{and}& B\left(\mathit{r}\right)+\omega \left(\mathit{r}\right)-\delta >0\end{array}& \begin{array}{cc}\text{outside}& S\end{array}\end{array}$$E6

Kalker [5] demonstrated that the solution to this contact problem can be found minimizing the total complementary energy $\left(V^{\ast }\right)$ 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^{\ast }=\frac{1}{2}{\int }_{S}p\left(\mathit{r}\right)\bullet \omega \left(\mathit{r}\right)\bullet \mathit{dS}+{\int }_{S}p\left(\mathit{r}\right)\bullet \left[B\left(\mathit{r}\right)-\delta \right]\bullet \mathit{dS}$$E7

To enable the numerical solution, the potential contact area is discretized into a set of $n$ pressure elements ${\mathrm{\Delta }}_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 ${\mathrm{\Delta }}_j$ is denoted by vector ${\mathit{r}}_{\mathit{j}}$.

Figure 4.

Discretization of the potential contact area into a pressure element mesh.

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

$$V^{\ast }=\frac{1}{2}\sum _{i=1}^n\left[p_i{\int }_{A_i}\omega \left(\mathit{r}\right)\bullet d{A}_i\right]+\sum _{i=1}^n\left[p_i{\int }_{A_i}\left[B\left(\mathit{r}\right)-\delta \right]\bullet d{A}_i\right]$$E8

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

$$\omega \left(\mathit{r}\right)=\sum _{j=1}^n\left[p_j\bullet f_j^{\left(1\right)}\left(\mathit{r}\right)\right]+\sum _{j=1}^n\left[p_j\bullet f_j^{\left(2\right)}\left(\mathit{r}\right)\right]=\sum _{j=1}^n\left[p_j\bullet t_j\left(\mathit{r}\right)\right]$$E9

where $t_j\left(\mathit{r}\right)$ is defined as the cumulated influence coefficient of pressure element ${\Delta }_j$ over point $Q$

$$t_j\left(\mathit{r}\right)=f_j^{\left(1\right)}\left(\mathit{r}\right)+f_j^{\left(2\right)}\left(\mathit{r}\right)$$E10

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

$$V^{\ast }=\frac{1}{2}\sum _{i=1}^n\left[p_i\bullet \sum _{j=1}^n{p}_j{\int }_{A_i}{t}_j\left(\mathit{r}\right)\bullet d{A}_i\right]+\sum _{i=1}^n\left[p_i{\int }_{A_i}\left[B\left(\mathit{r}\right)-\delta \right]\bullet d{A}_i\right]$$E11

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

1. The distance between the surfaces of the bodies $B\left(\mathit{r}\right)$ is assumed to be constant in all the pressure element ${\mathrm{\Delta }}_i$, and equal to the distance between the surface of the two bodies at its centroid, in such a way that $B\left(\mathit{r}\right)=B\left({\mathit{r}}_i\right)=B_i$.

2. The cumulated influence coefficient $t_j\left(\mathit{r}\right)$ is assumed to be constant over all the pressure element ${\mathrm{\Delta }}_i$, and equal to the value in its centroid, in such a way that $t_j\left(\mathit{r}\right)=t_j\left({\mathit{r}}_{\mathit{i}}\right)=t_{j,i}$. The coefficient $t_{j,i}$ can be defined as the cumulated influence coefficient of element ${\mathrm{\Delta }}_j$ over the centroid of element ${\mathrm{\Delta }}_i$, and it may be expressed as:

$$t_{j,i}=f_{j,i}^{\left(1\right)}+f_{j,i}^{\left(2\right)}$$

where $f_{j,i}^{\left(1\right)}$ and $f_{j,i}^{\left(2\right)}$ are the influence coefficients of element ${\Delta }_j$ over the centroid of element ${\Delta }_i$, which can be determined for each contacting body using Eq. (3).

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

$$V^{\ast }=\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{p}_i\bullet p_j\bullet t_{j,i}\bullet A_i+\sum _{i=1}^n{p}_i\bullet \left[B_i-\delta \right]\bullet A_i$$E12

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

$$\begin{array}{ccc}\frac{\partial V^{\ast }}{\partial p_i}=\sum _{j=1}^n{p}_j\bullet t_{j,i}\bullet A_i+\left[B_i-\delta \right]\bullet A_i=0& \mathit{\text{if}}& p_i>0\\ \frac{\partial V^{\ast }}{\partial p_i}=\sum _{j=1}^n{p}_j\bullet t_{j,i}\bullet A_i+\left[B_i-\delta \right]\bullet A_i\ge 0& \mathit{\text{if}}& p_i=0\end{array}$$E13

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 $\left(F_T\right)$ can be calculated as

$$F_T=\sum _{i=1}^n{p}_i\bullet A_i$$E14

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 $f_{j,i}$ for each body with finite dimensions in the following way:

$$f_{j,i}=f_{\mathit{jo},i}+\stackrel{\Psi }{\overbrace{\left[1.29-\frac{1}{1-\nu }\bullet \left(0.08-0.5\bullet \nu \right)\right]}}\bullet \sum _{m=1}^n{f}_{\mathit{jm},i}$$E15

where $f_{j0,i}$ is the influence coefficient of the original pressure element ${\Delta }_{\mathit{jo}}$, $\Psi$ is a correction factor proposed by Guilbault [7], and $f_{\mathit{jm},i}$ are the influence coefficients of the pressure elements ${\Delta }_{\mathit{jm}}$, that are mirrored instances of the element ${\Delta }_{\mathit{jo}}$ respect to the planes that coincide with the $n$ free-stress surfaces of the body.

Figure 5.

Contact between bodies of finite dimensions.

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 $\left(n+1\right)$, 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 $\left(\Gamma \right)$ 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).

Figure 6.

Example of a quadtree decomposition of the domain.

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 $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 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 $\left(\mathrm{\Gamma }\right)$ 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 $\Pi$, as shown in Figure 7.

Figure 7.

Definition of the potential contact area $\mathrm{\Gamma }.$

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 $\left(O_L\right)$ and a vector defining the contact normal.

3. The magnitude of the contact force $\left(F_T\right)$.

4. The initial level of uniform mesh density $\left(L_{\mathit{uni}}\right)$, which is a parameter that describes the size of the elements of the initial uniform pressure element mesh.

Figure 8.

Main algorithm of the proposed approach and algorithm to determine elements to split.

The algorithm starts determining the common tangent plane $\Pi$, where a local Cartesian coordinate system is defined, being the $Z_L$ axis normal to the plane $\Pi$ (step A1). The boundaries of the contacting bodies are normally projected onto the plane $\Pi$ to determine the potential contact area $\Gamma$ (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 $L_{\mathit{uni}}$ is reached for all the cells of the quadtree (step A3). All leaf cells of the quadtree are considered pressure elements ${\mathrm{\Delta }}_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 ${\mathrm{\Lambda }}_i=\mathit{\text{TRUE}}$, indicating that their properties (associated area $A_i$, normal gap $B_i$, contact pressure $p_i$, and cumulated influence coefficients $t_{j,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 $B_i$ for all the pressure elements (step A4) where ${\mathrm{\Lambda }}_i=\mathit{\text{TRUE}}$. The cumulated influence coefficients $t_{j,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 $p_i$ for each element ${\mathrm{\Delta }}_i$ of the discretization.

The flag ${\mathrm{\Lambda }}_i$ is defined as $\mathit{\text{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 ${\mathrm{\Lambda }}_i=\mathit{\text{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 $\left(N_f\right)$ can be determined afterwards using the following equation:

$$N_f=\sum _{i=1}^t\left[2\bullet n_{\left(i\right)}\bullet n_{\mathit{new}\left(i\right)}-n_{\mathit{new}\left(i\right)}^2\right]$$E16

where $t$ is the number of iterations performed by the algorithm, $n_{\left(i\right)}$ is the number of elements in iteration $i$, and $n_{\mathit{new}\left(i\right)}$ 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 $\lambda$) 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 ${\mathrm{\Delta }}_i$ will have an associated value ${\mathrm{\lambda }}_i$ of the observed physical magnitude, and in consequence, a discrete rate of change ${\phi }_{j,i}$ of $\lambda$ can be established between an element ${\mathrm{\Delta }}_i$ and any of its neighbors ${\mathrm{\Delta }}_j$ as

$${\phi }_{j,i}=\frac{\left|{\lambda }_j-{\lambda }_i\right|}{\mathit{max}\left(\left|{\lambda }_j\right|,\left|{\lambda }_i\right|\right)}$$E17

If the discrete rate of change ${\phi }_{j,i}$ between a pressure element ${\mathrm{\Delta }}_i$ and any of its neighbors ${\mathrm{\Delta }}_j$ is higher than an arbitrarily defined value ${\phi }_{\mathit{max}}$ (representing the maximum allowed rate of change of the observed physical magnitude), then the pressure element ${\mathrm{\Delta }}_i$ is marked as a candidate to be split.

However, in some situations, the rate of change of $\lambda$ is so high that the condition ${\phi }_{j,i}<{\phi }_{\mathit{max}}$ cannot be reached, and the refinement strategy based on the rate of change of $\lambda$ 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 $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_{\mathit{max}}$, referred to as the maximum degree of mesh refinement.

It is important to point out that ${\phi }_{\mathit{max}}$ is a target value and may not be always reached. If the maximum degree of mesh refinement $L_{\mathit{max}}$ is reached for all pressure elements before the target value for ${\phi }_{\mathit{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 $A_i$, normal gap $B_i$, contact pressure $p_i$, and cumulated influence coefficients $t_{j,i}$).

2. The quadtree data structure.

3. The maximum degree of mesh refinement $\left(L_{\mathit{max}}\right)$.

4. The maximum allowed rate of change of the observed physical magnitude $\left({\phi }_{\mathit{max}}\right)$.

The algorithm starts defining the flag $K_i$ as $\mathit{\text{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 ${\mathrm{\Delta }}_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 ${\mathrm{\Delta }}_i$ and ${\mathrm{\Delta }}_j$, the algorithm determines the associated physical magnitudes ${\lambda }_i$ and ${\lambda }_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 ${\phi }_{j,i}$ of the observed magnitude between pressure element ${\mathrm{\Delta }}_i$ and his neighbor ${\mathrm{\Delta }}_j$ is obtained using Eq. (17) (step B3). If ${\phi }_{j,i}$ is lower than the user-defined value ${\phi }_{\mathit{max}}$, the next neighbor pressure element ${\mathrm{\Delta }}_{j+1}$ is evaluated. In contrast, if the discrete rate of change of the observed magnitude between both elements is greater than ${\phi }_{\mathit{max}}$, then element ${\mathrm{\Delta }}_i$ is considered as a candidate to be split. Two additional conditions must be fulfilled so ${\mathrm{\Delta }}_i$ can be marked to be split:

1. On one hand, $L_i$ must be lower than $L_{\mathit{max}}$, to avoid that the algorithm refines the mesh indefinitely in those cases where ${\phi }_{\mathit{max}}$ cannot be reached.

2. On the other hand, $L_i$ must be less than or equal to $L_j$, 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 ${\mathrm{\Delta }}_i$ is marked to be split by defining ${\mathrm{K}}_i=\mathit{\text{TRUE}}$. The algorithm finishes when all the pressure elements have been evaluated, returning an array that contains the indices of those elements where ${\mathrm{K}}_i=\mathit{\text{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 (${\lambda }_i=p_i$), 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 ${\mathrm{\Delta }}_i$ that is within the contact area ($p_i>0$) and of an adjacent element ${\mathrm{\Delta }}_j$ that is outside of the contact area ($p_j=0$) is always ${\phi }_{j,i}=1$. Therefore, if a value lower than 1 is specified for ${\phi }_{\mathit{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, $L_{\mathit{uni}}$.

2. The maximum degree of mesh refinement, $L_{\mathit{max}}$.

3. The maximum allowed rate of change of the physical magnitude, ${\phi }_{\mathit{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 ($L_{\mathit{uni}}\ge L_{\mathit{max}}$, ${\phi }_{\mathit{max}}$ not relevant): using this setting, the contact problem is solved using a uniform mesh, whose mesh density is defined by $L_{\mathit{uni}}$.

2. Setting 2 ($L_{\mathit{uni}}<L_{\mathit{max}},{\phi }_{\mathit{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 $L_{\mathit{max}}$.

3. Setting 3 ($L_{\mathit{uni}}<L_{\mathit{max}},{\phi }_{\mathit{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.

Figure 9.

Definition of the indenters for (a) case of study I and (b) case of study II.

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

In both cases, the root cell of the quadtree results in a $20\times 20\mathrm{mm}$ 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 $L_{\mathit{uni}}$ has been varied, keeping $L_{\mathit{uni}}=L_{\mathit{max}}$ and ${\phi }_{\mathit{max}}=0$ (setting 1 in Section 3.1). Figure 10a–c 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.

Figure 10.

Axisymmetric representation of the resulting contact area and pressure element mesh obtained for CoSI under several configurations of the approach.

The contact pressure distributions along the principal axes of the contact area of the solutions shown in Figure 10a–c 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 $L_{\mathit{uni}}$), the results obtained by the proposed approach converge toward the reference solution.

Figure 11.

Contact pressure distribution for CoSI under several configurations of the approach.

Using this configuration of the approach, a mesh containing $4^{L_{\mathit{uni}}}$ pressure elements is used, regardless of the nature of the contact problem to be solved. Under these circumstances, the computational cost is proportional to $4^{2\bullet L_{\mathit{uni}}}$, 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 $L_{\mathit{uni}}$, 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 $L_{\mathit{uni}}$ and $L_{\mathit{max}}$ have been varied, keeping $L_{\mathit{uni}}<L_{\mathit{max}}$ and ${\phi }_{\mathit{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.

Figure 12.

Axisymmetric representation of the resulting contact area and pressure element mesh obtained for CoSII under two different configurations 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 $L_{\mathit{uni}}$. For any given value of $L_{\mathit{max}}$, 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 $L_{\mathit{uni}}$. 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 $L_{\mathit{max}}$ and $L_{\mathit{uni}}$. 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 ${\phi }_{\mathit{max}}$ has been varied, keeping $L_{\mathit{uni}}<L_{\mathit{max}}$ (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 10d–f 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 ${\phi }_{\mathit{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.

Figure 13.

Contact pressure distribution for CoSII under several configurations of the approach.

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 ${\phi }_{\mathit{max}}$, since the same values are obtained regardless of the value selected for this parameter. This is because when ${\phi }_{\mathit{max}}<1$, the accuracy in which the border of the contact area is computed depends only on the value selected for $L_{\mathit{max}}$, 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 ${\phi }_{\mathit{max}}>0$. Although this reduction is not as important as the one achieved by maximizing $L_{\mathit{max}}-L_{\mathit{uni}}$ (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 $L_{\mathit{uni}}$), 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 ${\phi }_{\mathit{max}}$ and (ii) the maximum degree of mesh refinement, defined by the parameter $L_{\mathit{max}}$.

The configuration of the approach is defined by a unique combination of values for $L_{\mathit{uni}}$, $L_{\mathit{max}}$, and ${\phi }_{\mathit{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 $L_{\mathit{uni}}=L_{\mathit{max}}$, a uniform mesh is used to solve the contact problem, regardless of the value selected for ${\phi }_{\mathit{max}}$. Under this configuration, it can be observed that the obtained results converge toward the reference solution as $L_{\mathit{uni}}$ is increased. However, an exponential growth of the computational cost is produced as the pressure element mesh is refined.

2. When $L_{\mathit{uni}}<L_{\mathit{max}}$ and ${\phi }_{\mathit{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 $L_{\mathit{max}}-L_{\mathit{ini}}$, while the accuracy of the solution remains unaffected.

3. In last place, when $L_{\mathit{uni}}<L_{\mathit{max}}$ and ${\phi }_{\mathit{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 ${\phi }_{\mathit{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]).