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:


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

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

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

69

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

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 Δ<sup>i</sup> for the initial iteration of the algorithm, which is performed using a uniform pressure element mesh. All the elements are marked with the flag Λ<sup>i</sup> ¼ TRUE, indicating that their properties

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

, and cumulated influence coefficients tj,i)

potential contact area Γ (step A2).

are not computed yet.

(associated area Ai, normal gap Bi, contact pressure pi

Figure 7. Definition of the potential contact area Γ:

Adaptive Mesh Refinement Strategy for the Semi-Analytical Solution of Frictionless Elastic Contact Problems http://dx.doi.org/10.5772/intechopen.72422 69

3. Computational approach to solve frictionless elastic contact problems

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

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 algo-

The main algorithm of the approach to solve contact problems using adaptive mesh refinement

iv. The initial level of uniform mesh density ð Þ Luni , which is a parameter that describes the

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

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

intersection of the projection of the bodies on the plane Π, as shown in Figure 7.

rithms are computationally efficient and easy to implement.

iii. The magnitude of the contact force ð Þ FT .

Figure 7. Definition of the potential contact area Γ:

is shown in Figure 8. The following inputs are required by the algorithm:

ii. The initial point of contact ð Þ OL and a vector defining the contact normal.

size of the elements of the initial uniform pressure element mesh.

using adaptive mesh refinement

68 Contact and Fracture Mechanics

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

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 Δ<sup>i</sup> for the initial iteration of the algorithm, which is performed using a uniform pressure element mesh. All the elements are marked with the flag Λ<sup>i</sup> ¼ 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 Λ<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 pi for each element Δ<sup>i</sup> of the discretization.

The flag Λ<sup>i</sup> 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 Λ<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 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:

$$N\_f = \sum\_{i=1}^{t} \left[ 2 \cdot n\_{(i)} \cdot n\_{new(i)} - n\_{new(i)}^2 \right] \tag{16}$$

change of the observed physical magnitude), then the pressure element Δ<sup>i</sup> is marked as a

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

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

71

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

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

The main routine of the algorithm to determine elements to split is shown in Figure 8b. The

i. The properties associated with the pressure elements (area Ai, normal gap Bi, contact

iv. The maximum allowed rate of change of the observed physical magnitude wmax ð Þ.

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

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

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

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

, and cumulated influence coefficients tj,i).

candidate to be split.

pressure pi

ii. The quadtree data structure.

maximum degree of mesh refinement.

value for wmax is achieved, the mesh refinement will finish.

following input information is required by the algorithm:

iii. The maximum degree of mesh refinement ð Þ Lmax .

assumed that none of the elements will be divided.

the indexes of the k neighbors of a given pressure element.

determined from these values by performing additional calculations.

indefinitely in those cases where wmax cannot be reached.

where t is the number of iterations performed by the algorithm, nð Þ<sup>i</sup> is the number of elements in iteration i, and nnew ið Þ is the number of new elements in iteration i.
