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

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

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,

ωð Þ¼ r f <sup>j</sup>

<sup>þ</sup> <sup>A</sup>∙ln <sup>C</sup> <sup>þ</sup>

A ¼ dy þ b C ¼ dx þ a B ¼ dy � b D ¼ dx � a

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.

ð Þr is the influence coefficient of pressure element Δ<sup>j</sup> over the point H, which can be

D þ

. The area of the pressure element shown in Figure 1b is Aj ¼ 2a � 2b,

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>A</sup><sup>2</sup> <sup>þ</sup> <sup>C</sup><sup>2</sup> <sup>p</sup>

" #

" #) (3)

( " #

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>A</sup><sup>2</sup> <sup>þ</sup> <sup>D</sup><sup>2</sup> <sup>p</sup> . Under these condi-

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>B</sup><sup>2</sup> <sup>þ</sup> <sup>D</sup><sup>2</sup> <sup>p</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>A</sup><sup>2</sup> <sup>þ</sup> <sup>D</sup><sup>2</sup> <sup>p</sup>

ð Þr ∙pj (2)

<sup>þ</sup> <sup>D</sup>∙ln <sup>B</sup> <sup>þ</sup>

A þ

with a specific shape (such as triangles, rectangles, hexagons, etc.).

and the uniform pressure distribution that acts over this area is p r<sup>0</sup> ð Þ¼ pj

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>A</sup><sup>2</sup> <sup>þ</sup> <sup>C</sup><sup>2</sup> <sup>p</sup>

" #

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>B</sup><sup>2</sup> <sup>þ</sup> <sup>D</sup><sup>2</sup> <sup>p</sup>

and will be denoted by Δ<sup>j</sup>

62 Contact and Fracture Mechanics

analytically determined as

1 � ν

where f <sup>j</sup>

f j ð Þ¼ r

tions, the closed-form solution for Eq. (1) is

<sup>2</sup>π<sup>G</sup> <sup>C</sup>∙ln <sup>A</sup> <sup>þ</sup>

C þ

<sup>þ</sup>B∙ln <sup>D</sup> <sup>þ</sup>

B þ

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

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>B</sup><sup>2</sup> <sup>þ</sup> <sup>D</sup><sup>2</sup> <sup>p</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>B</sup><sup>2</sup> <sup>þ</sup> <sup>C</sup><sup>2</sup> <sup>p</sup>

Figure 2. Normal deflection produced by a complex pressure distribution.

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 r xL; yL; zL � �, being zL <sup>¼</sup> 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 pð Þr is generated in the true contact area S that deforms the contacting bodies. In this way, elastic normal deflections are produced

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

in the surfaces of the bodies 1 and 2, which are denoted by <sup>ω</sup>ð Þ<sup>1</sup> ð Þ<sup>r</sup> and <sup>ω</sup>ð Þ<sup>2</sup> ð Þ<sup>r</sup> , respectively. The total normal deflection is denoted by ωð Þr and can be calculated as

$$
\omega(\mathbf{r}) = \omega^{(1)}(\mathbf{r}) + \omega^{(2)}(\mathbf{r}) \tag{5}
$$

<sup>ω</sup>ð Þ¼ <sup>r</sup> <sup>X</sup><sup>n</sup>

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

<sup>V</sup><sup>∗</sup> <sup>¼</sup> <sup>1</sup> 2 Xn i¼1

where f

ð Þ1 j,i and f

calculated as

ð Þ2

j¼1 pj ∙f ð Þ1 <sup>j</sup> ð Þr h i þX<sup>n</sup>

pi ∙ Xn j¼1 pj ð Ai

2 4

centroid, in such a way that Bð Þ¼ r Bð Þ¼ r<sup>i</sup> Bi.

over the centroid of element Δi, and it may be expressed as:

which can be determined for each contacting body using Eq. (3).

<sup>V</sup><sup>∗</sup> <sup>¼</sup> <sup>1</sup> 2 Xn i¼1

minimizing Eq. (12) under the following restrictions:

∂V<sup>∗</sup> ∂pi

> ∂V<sup>∗</sup> ∂pi

<sup>¼</sup> <sup>X</sup><sup>n</sup> j¼1 pj

<sup>¼</sup> <sup>X</sup><sup>n</sup> j¼1 pj j¼1 pj ∙f ð Þ2 <sup>j</sup> ð Þr h i <sup>¼</sup> <sup>X</sup><sup>n</sup>

where tjð Þr is defined as the cumulated influence coefficient of pressure element Δ<sup>j</sup> over point Q

ð Þ2

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

pi ð Ai

ð Þ1 <sup>j</sup> ð Þþ r f

tjð Þr ∙dAi

i. The distance between the surfaces of the bodies Bð Þr 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

ii. The cumulated influence coefficient tjð Þr is assumed to be constant over all the pressure element Δi, and equal to the value in its centroid, in such a way that tjð Þ¼ r tjð Þ¼ ri tj,i. The coefficient tj,i can be defined as the cumulated influence coefficient of element Δ<sup>j</sup>

j,i are the influence coefficients of element Δ<sup>j</sup> over the centroid of element Δi,

i¼1 pi

∙tj,i∙Ai þ ½ � Bi � δ ∙Ai ¼ 0 if pi > 0

∙tj,i∙Ai þ ½ � Bi � δ ∙Ai ≥ 0 if pi ¼ 0

<sup>∙</sup>tj,i∙Ai <sup>þ</sup>X<sup>n</sup>

The solution to the contact problem, in terms of contact pressure distribution, can be found by

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

tj,i ¼ f ð Þ1 j,i þ f ð Þ2 j,i

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

Xn j¼1 pi ∙pj 3 <sup>5</sup> <sup>þ</sup>X<sup>n</sup> i¼1

tjð Þ¼ r f

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

j¼1 pj ∙tjð Þr

h i (9)

65

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

<sup>j</sup> ð Þr (10)

� � (11)

∙½ � Bi � δ ∙Ai (12)

(13)

½ � Bð Þ� r δ ∙dAi

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{cases} p(\mathbf{r}) > 0 \text{ and } \; B(\mathbf{r}) + \omega(\mathbf{r}) - \delta = 0 & \text{inside} \quad \mathbf{S} \\ p(\mathbf{r}) = 0 \quad \text{and} \quad B(\mathbf{r}) + \omega(\mathbf{r}) - \delta > 0 & \text{outside} \quad \mathbf{S} \end{cases} \tag{6}$$

Kalker [5] demonstrated that the solution to this contact problem can be found minimizing the total complementary energy <sup>V</sup><sup>∗</sup> ð Þ 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^\* = \frac{1}{2} \int\_S p(\mathbf{r}) \cdot \omega(\mathbf{r}) \cdot d\mathbf{S} + \int\_S p(\mathbf{r}) \cdot [B(\mathbf{r}) - \delta] \cdot d\mathbf{S} \tag{7}$$

To enable the numerical solution, the potential contact area is discretized into a set of n pressure elements Δ<sup>j</sup> (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 Δ<sup>j</sup> is denoted by vector rj.

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

$$V^\* = \frac{1}{2} \sum\_{i=1}^n \left[ p\_i \int\_{A\_i} \omega(\mathbf{r}) \cdot dA\_i \right] + \sum\_{i=1}^n \left[ p\_i \int\_{A\_i} [B(\mathbf{r}) - \delta] \cdot dA\_i \right] \tag{8}$$

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

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

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

$$\omega(\mathbf{r}) = \sum\_{j=1}^{n} \left[ p\_j \mathbf{\dot{f}}\_j^{(1)}(\mathbf{r}) \right] + \sum\_{j=1}^{n} \left[ p\_j \mathbf{\dot{f}}\_j^{(2)}(\mathbf{r}) \right] = \sum\_{j=1}^{n} \left[ p\_j \mathbf{\dot{t}}\_j(\mathbf{r}) \right] \tag{9}$$

where tjð Þr is defined as the cumulated influence coefficient of pressure element Δ<sup>j</sup> over point Q

$$t\_{\dot{j}}(\mathbf{r}) = f\_{\dot{j}}^{(1)}(\mathbf{r}) + f\_{\dot{j}}^{(2)}(\mathbf{r}) \tag{10}$$

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

in the surfaces of the bodies 1 and 2, which are denoted by <sup>ω</sup>ð Þ<sup>1</sup> ð Þ<sup>r</sup> and <sup>ω</sup>ð Þ<sup>2</sup> ð Þ<sup>r</sup> , respectively. The

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

pð Þr > 0 and Bð Þþ r ωð Þ� r δ ¼ 0 inside S

Kalker [5] demonstrated that the solution to this contact problem can be found minimizing the total complementary energy <sup>V</sup><sup>∗</sup> ð Þ 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

> ð S

þX<sup>n</sup> i¼1

pi ð Ai

½ � Bð Þ� r δ ∙dAi � �

To enable the numerical solution, the potential contact area is discretized into a set of n pressure elements Δ<sup>j</sup> (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

pð Þr ∙ωð Þr ∙dS þ

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

ωð Þr ∙dAi � �

pi ð Ai

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

<sup>ω</sup>ð Þ¼ <sup>r</sup> <sup>ω</sup>ð Þ<sup>1</sup> ð Þþ <sup>r</sup> <sup>ω</sup>ð Þ<sup>2</sup> ð Þ<sup>r</sup> (5)

pð Þr ∙½ � Bð Þ� r δ ∙dS (7)

(8)

<sup>p</sup>ð Þ¼ <sup>r</sup> 0 and <sup>B</sup>ð Þþ <sup>r</sup> <sup>ω</sup>ð Þ� <sup>r</sup> <sup>δ</sup> <sup>&</sup>gt; <sup>0</sup> outside <sup>S</sup> (6)

total normal deflection is denoted by ωð Þr and can be calculated as

<sup>V</sup><sup>∗</sup> <sup>¼</sup> <sup>1</sup> 2 ð S

pressure element Δ<sup>j</sup> is denoted by vector rj.

<sup>V</sup><sup>∗</sup> <sup>¼</sup> <sup>1</sup> 2 Xn i¼1

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

size are unknown

64 Contact and Fracture Mechanics

complementary energy as

$$V^\* = \frac{1}{2} \sum\_{i=1}^n \left[ p\_i \cdot \sum\_{j=1}^n p\_j \int\_{A\_i} t\_j(\mathbf{r}) \cdot dA\_i \right] + \sum\_{i=1}^n \left[ p\_i \int\_{A\_i} [B(\mathbf{r}) - \delta] \cdot dA\_i \right] \tag{11}$$

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


$$t\_{\dot{\mu}^i} = f\_{\dot{\mu}^i}^{(1)} + f\_{\dot{\mu}^i}^{(2)}$$

where f ð Þ1 j,i and f ð Þ2 j,i are the influence coefficients of element Δ<sup>j</sup> 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^\* = \frac{1}{2} \sum\_{i=1}^n \sum\_{j=1}^n p\_i \cdot p\_j \cdot t\_{\flat, i} \cdot A\_i + \sum\_{i=1}^n p\_i \cdot [B\_i - \delta] \cdot A\_i \tag{12}$$

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

$$\begin{aligned} \frac{\partial V^\*}{\partial p\_i} &= \sum\_{j=1}^n p\_j \cdot t\_{j,i} \cdot A\_i + [B\_i - \delta] \cdot A\_i = 0 \quad \text{if} \quad p\_i > 0\\ \frac{\partial V^\*}{\partial p\_i} &= \sum\_{j=1}^n p\_j \cdot t\_{j,i} \cdot A\_i + [B\_i - \delta] \cdot A\_i \ge 0 \quad \text{if} \quad p\_i = 0 \end{aligned} \tag{13}$$

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

$$F\_T = \sum\_{i=1}^{n} p\_i \cdot A\_i \tag{14}$$

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

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

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

Figure 5. Contact between bodies of finite dimensions.

Figure 6. Example of a quadtree decomposition of the domain.

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.

, and

67

#### 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\_{j\alpha,i} + \overbrace{\left[1.29 - \frac{1}{1-\nu} \cdot (0.08 - 0.5 \cdot \nu)\right]}^{\Psi} \cdot \sum\_{m=1}^{n} f\_{jm,i} \tag{15}$$

where f <sup>j</sup>0,i is the influence coefficient of the original pressure element Δjo, Ψ is a correction factor proposed by Guilbault [7], and f jm,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).

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

Figure 5. Contact between bodies of finite dimensions.

FT <sup>¼</sup> <sup>X</sup><sup>n</sup> i¼1 pi

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

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

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

1 � ν

� � zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ Ψ

where f <sup>j</sup>0,i is the influence coefficient of the original pressure element Δjo, Ψ is a correction factor proposed by Guilbault [7], and f jm,i are the influence coefficients of the pressure elements Δjm, that are mirrored instances of the element Δjo respect to the planes that coincide

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

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

∙ð Þ 0:08 � 0:5∙ν

∙ Xn m¼1

f jm,i (15)

the described method can lead to erroneous solutions of the contact problem.

method to consider contact bodies of finite dimensions is described in this section.

<sup>f</sup> j,i <sup>¼</sup> <sup>f</sup> jo,i <sup>þ</sup> <sup>1</sup>:<sup>29</sup> � <sup>1</sup>

the number of finite dimensions taken into account in the problem.

with the n free-stress surfaces of the body.

2.4. Quadtree decomposition of the domain

2.3. Contact between bodies of finite dimensions

the following way:

66 Contact and Fracture Mechanics

South-East (SE).

∙Ai (14)

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.

Figure 6. Example of a quadtree decomposition of the domain.
