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 Δ<sup>i</sup> will have an associated value λ<sup>i</sup> of the observed physical magnitude, and in consequence, a discrete rate of change wj,i of λ can be established between an element Δ<sup>i</sup> and any of its neighbors Δ<sup>j</sup> as

$$\varphi\_{j,i} = \frac{|\lambda\_j - \lambda\_i|}{\max\left(|\lambda\_j|, |\lambda\_i|\right)}\tag{17}$$

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

However, in some situations, the rate of change of λ is so high that the condition wj,i < wmax 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 wmax 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 wmax 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:


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

Then, an iterative process starts whose first step is the determination of the normal gap Bi for all the pressure elements (step A4) where Λ<sup>i</sup> ¼ 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

The flag Λ<sup>i</sup> is defined as FALSE for all the elements present in the discretization, indicating that

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 Λ<sup>i</sup> ¼ 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

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

� � can be determined afterwards using the following equation:

new ið Þ

�; <sup>λ</sup><sup>i</sup> j j � � (17)

(16)

<sup>2</sup>∙nð Þ<sup>i</sup> <sup>∙</sup>nnew ið Þ � <sup>n</sup><sup>2</sup>

where t is the number of iterations performed by the algorithm, nð Þ<sup>i</sup> is the number of elements

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 Δ<sup>i</sup> will have an associated value λ<sup>i</sup> of the observed physical magnitude, and in consequence, a discrete rate of change wj,i of λ can be established between an element Δ<sup>i</sup> and any of its neighbors Δ<sup>j</sup> as

> <sup>w</sup>j,i <sup>¼</sup> <sup>λ</sup><sup>j</sup> � <sup>λ</sup><sup>i</sup> � � � �

max λ<sup>j</sup> � � �

If the discrete rate of change wj,i between a pressure element Δ<sup>i</sup> and any of its neighbors Δ<sup>j</sup> is higher than an arbitrarily defined value wmax (representing the maximum allowed rate of

h i

every point in the region [9].

70 Contact and Fracture Mechanics

proposed approach Nf

pi for each element Δ<sup>i</sup> of the discretization.

the properties of these elements have already been computed.

starts again (step A4), and it is repeated until no new elements are created.

Nf <sup>¼</sup> <sup>X</sup><sup>t</sup> i¼1

in iteration i, and nnew ið Þ is the number of new elements in iteration i.

3.1. Algorithm to determine elements to split


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 Δ<sup>i</sup> 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 Δ<sup>i</sup> and Δj, the algorithm determines the associated physical magnitudes λ<sup>i</sup> and λ<sup>j</sup> (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 wj,i of the observed magnitude between pressure element Δ<sup>i</sup> and his neighbor Δ<sup>j</sup> is obtained using Eq. (17) (step B3). If wj,i is lower than the user-defined value wmax, the next neighbor pressure element Δ<sup>j</sup>þ<sup>1</sup> is evaluated. In contrast, if the discrete rate of change of the observed magnitude between both elements is greater than wmax, then element Δ<sup>i</sup> is considered as a candidate to be split. Two additional conditions must be fulfilled so Δ<sup>i</sup> can be marked to be split:

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

ii. 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.

i. Setting 1 (Luni ≥ Lmax, wmax not relevant): using this setting, the contact problem is solved

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

http://dx.doi.org/10.5772/intechopen.72422

73

ii. Setting 2 (Luni < Lmax, wmax ¼ 0): using this setting, the contact problem is solved using adaptive mesh refinement outside the true contact area. Inside the true contact area, a

iii. Setting 3 (Luni < Lmax, wmax > 0): using this setting, the contact problem is solved using

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

The performance of the proposed approach is illustrated in this section, considering its accu-

• 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

• 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

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

In both cases, the root cell of the quadtree results in a 20 � 20 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

racy and computational cost. For such a purpose, two cases of study are considered:

mechanical components such as ball bearings, gears and rail-wheel systems.

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.

adaptive mesh refinement both inside and outside the true contact area.

using a uniform mesh, whose mesh density is defined by Luni.

uniform mesh is used, whose mesh density is defined by Lmax.

Eq. (16).

considered.

4. Numerical examples

method described in Section 2.3.

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