Open access peer-reviewed chapter

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

By Victor Roda-Casanova and Francisco Sanchez-Marin

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

DOI: 10.5772/intechopen.72422

Downloaded: 214

Abstract

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

Keywords

  • contact analysis
  • semi-analytical methods
  • adaptive refinement

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 npressure elements, with a uniform pressure distribution assumed to be acting over each one of them. Influence coefficients are used to relate the pressure applied over each pressure element with the displacements that this pressure produces at the centroid of the other elements of the mesh. Using these influence coefficients, the solution to the contact problem can be numerically found in terms of the contact pressure distribution that satisfies the contact conditions.

As usual, in numerical methods based on the discretization of the domain, the election of the number of pressure elements in which the domain is divided involves a commitment between accuracy and computational cost. Kalker [1] stated that the computational cost of these semi-analytical methods can be defined by the number of influence coefficients that need to be calculated to solve the contact problem (that, in general, is proportional to n2). He also argued that the accuracy of the solution to the contact problem, in terms of contact area and contact pressure distribution, depends on the refinement of the pressure element mesh, especially in those regions close to the border of the contact area. Consequently, an improvement of the accuracy of the results necessarily implies an increment of the computational cost.

When both shape and location of the true contact area are known in advance, the efficiency of the method can be maximized by discretizing an area similar to the true contact area. But when the true contact area is unknown, it is difficult to optimize the efficiency of the method, since the whole potential contact area must be discretized to consider any possible shape and location of the true contact area. In those cases, it is common to use a uniform pressure element mesh for the whole domain, being more or less dense depending on the desired accuracy and on the capabilities of the computer used to solve the contact problem. In consequence, there could be many pressure elements in the discretization out of the true contact area, what causes a loss in the efficiency of the method.

These difficulties could be partially overcome using adaptive mesh refinement strategies. These techniques have been previously used to improve the efficiency of numerical methods based in the discretization of the domain, especially in FEM procedures [2]. However, no previous use of adaptive refinement has been found in the literature for the solution of contact problems using semi-analytical methods.

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

2. Theoretical background

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

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

Consider a body that, because of its main features, can be approached to an elastic half-space, as the one shown in Figure 1a. A Cartesian coordinate system is defined over the surface of this body, which Xand Yaxes define a plane that is coincident with its surface, and the Zaxis points inward him. A normal pressure distribution pis 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 Cwithin the area S, whose position is defined by the vector rxyz, being z=0. Consider another point Hin the surface of the body, whose position is defined by the vector rxyz, being z=0. The normal elastic deflection produced at a point Hdue to a normal pressure distribution applied over the area Scan be determined by the superposition of the Boussinesq relation [3]:

ωr=1ν2πGSprrrdSE1

where νis the Poisson coefficient and Gis 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 Sand on the considered pressure distribution. However, several closed-form solutions can be found in the literature for certain pressure distributions applied over areas with a specific shape (such as triangles, rectangles, hexagons, etc.).

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

ωr=fjrpjE2

where fjris the influence coefficient of pressure element jover the point H, which can be analytically determined as

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

where coefficients A, B, C, and Dare calculated as

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

These pressure elements can be useful to determine the normal displacement produced at the surface of a body due to a non-uniform pressure distribution applied over a complex area. To illustrate this methodology, consider a complex area S, as the one shown in Figure 2a, over which an arbitrary pressure distribution is acting. To determine the displacement field produced by this pressure distribution, the area Sis discretized into a mesh of nrectangular pressure elements Δj, as shown in Figure 2b. Then, the arbitrary pressure distribution is approached by assigning a uniform pressure value pjto 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

ωr=j=1npjfjrE4

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

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

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

Consider two bodies 1and 2in 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.

Figure 3.

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

A Cartesian coordinate system is defined with origin at point OL, being the local axis ZLnormal to the plane Πand pointing inward the body 2. Consider a generic point Qin the plane Π, whose position is defined by the vector rxLyLzL, being zL=0. The gap between the two bodies, measured along ZLaxis, 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 pris generated in the true contact area Sthat deforms the contacting bodies. In this way, elastic normal deflections are produced in the surfaces of the bodies 1and 2, which are denoted by ω1rand ω2r, respectively. The total normal deflection is denoted by ωrand can be calculated as

ωr=ω1r+ω2rE5

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

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

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

V=12SprωrdS+SprBrδdSE7

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

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=12i=1npiAiωrdAi+i=1npiAiBrδdAiE8

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

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

where tjris defined as the cumulated influence coefficient of pressure element Δjover point Q

tjr=fj1r+fj2rE10

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

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

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

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

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

tj,i=fj,i1+fj,i2

where fj,i1and fj,i2are the influence coefficients of element Δjover the centroid of element Δi, which can be determined for each contacting body using Eq. (3).

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

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

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

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

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

FT=i=1npiAiE14

2.3. Contact between bodies of finite dimensions

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

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

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

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

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

where fj0,iis the influence coefficient of the original pressure element Δjo, Ψis a correction factor proposed by Guilbault [7], and fjm,iare the influence coefficients of the pressure elements Δjm, that are mirrored instances of the element Δjorespect to the planes that coincide with the nfree-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 n+1, being nthe 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).

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 jin the structure is denoted by Lj, and represents the number of divisions performed from the root of the quadtree. According to this definition, Ljis 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 Ljis known. The root cell of the quadtree is usually denoted by level 0. Any cell that is not subdivided anymore is a leaf cell (displayed in gray in Figure 6e), while subdivided cells are referred to as non-leaf cells.

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

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

Figure 7.

Definition of the potential contact area Γ.

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 OLand a vector defining the contact normal.

  3. The magnitude of the contact force FT.

  4. The initial level of uniform mesh density Luni, which is a parameter that describes the size of the elements of the initial uniform pressure element mesh.

Figure 8.

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

The algorithm starts determining the common tangent plane Π, where a local Cartesian coordinate system is defined, being the ZLaxis normal to the plane Π(step A1). The boundaries of the contacting bodies are normally projected onto the plane Πto determine the potential contact area Γ(step A2).

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

To maximize the efficiency of the proposed approach, it is important to minimize the number of pressure elements located outside of the potential contact area. This can be achieved by ensuring that the potential contact area is enclosed by a root cell of the quadtree coincident with the minimum bounding rectangle (MBR), defined as the minimum rectangle that contains every point in the region [9].

Then, an iterative process starts whose first step is the determination of the normal gap Bifor all the pressure elements (step A4) where Λi=TRUE. The cumulated influence coefficients tj,iof 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 pifor each element Δiof the discretization.

The flag Λiis defined as FALSEfor all the elements present in the discretization, indicating that the properties of these elements have already been computed.

At this point, the adaptive mesh refinement is performed. For such a purpose, the algorithm that determines the elements that must be split (that is described in section 3.1) is called (step A7), which returns an array with the indexes of these elements. Then, the selected elements are split (step A8) and the quadtree data structure is updated with the information of the new elements, which are marked with the flag Λi=TRUE, indicating that their properties are not computed yet. If no new elements are created, the iterative process finishes and the contact results are displayed (step A9). In contrast, if new elements are created, the iterative process starts again (step A4), and it is repeated until no new elements are created.

The main advantage of this implementation is that only the normal gap and the influence coefficients related to the new elements are computed for each iteration, decreasing the global computational cost of the method. The number of influence coefficients calculated in the proposed approach Nfcan be determined afterwards using the following equation:

Nf=i=1t2ninnewinnewi2E16

where tis the number of iterations performed by the algorithm, niis the number of elements in iteration i, and nnewiis the number of new elements in iteration i.

3.1. Algorithm to determine elements to split

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

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

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

However, in some situations, the rate of change of λis so high that the condition φj,i<φmaxcannot 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 Ljthat 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 φmaxis a target value and may not be always reached. If the maximum degree of mesh refinement Lmaxis reached for all pressure elements before the target value for φmaxis achieved, the mesh refinement will finish.

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

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

  2. The quadtree data structure.

  3. The maximum degree of mesh refinement Lmax.

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

The algorithm starts defining the flag Kias FALSEfor all the elements present in the current discretization. The flag Kiindicates 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 kneighbors of every pressure element Δiin 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 kneighbors of a given pressure element.

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

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

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

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

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

3.2. Final remarks

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

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

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

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

The topology of the resulting pressure element mesh, inside and outside the true contact area, depends on the configuration of the proposed approach, which is defined by a unique combination of the three input parameters:

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

  2. The maximum degree of mesh refinement, Lmax.

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

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

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

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

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

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

4. Numerical examples

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

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

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

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 70GPaand a Poisson coefficient of 0.35. A total contact load FT=60kNis considered.

In both cases, the root cell of the quadtree results in a 20×20mmsquare. 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 Lunihas been varied, keeping Luni=Lmaxand φmax=0(setting 1 in Section 3.1). Figure 10ac show examples of the resulting contact area and pressure element mesh that have been obtained for case of study I under this setting of the approach. The computational cost of the proposed approach to solve the case of study I is also shown for each configuration.

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 10ac are shown in Figure 11a. As expected, it can be observed that as the pressure element mesh is refined (by increasing the value selected for Luni), the results obtained by the proposed approach converge toward the reference solution.

Figure 11.

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

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

4.2. Performance of the approach when adaptive mesh refinement is performed outside the true contact area

In this section, the performance of the proposed approach when adaptive refinement is performed outside the true contact area is illustrated. To do so, the contact problems defined by cases of study I and II are solved under several configurations of the approach, in which Luniand Lmaxhave been varied, keeping Luni<Lmaxand φ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 Luni. For any given value of Lmax, the same contact pressure distributions as the ones obtained with a uniform pressure element mesh for the whole domain have been obtained (shown in Figure 11a), regardless of the value selected for Luni. This implies that the variation of the pressure element mesh outside the true contact area does not have any impact on the solution of the contact problem.

On the other hand, comparing the computational cost of the solutions shown in Figure 10c (uniform mesh) and 10d (adaptive refinement outside the true contact area), it can be observed that an important reduction of the computational cost is achieved by increasing the difference between Lmaxand Luni. Similar tendencies are observed for case of study II, where the reductions of computational cost are even more remarkable due to the presence of finite dimensions.

4.3. Performance of the approach when adaptive mesh refinement is performed both inside and outside the true contact area

In this section, the performance of the proposed approach when adaptive refinement is also performed inside the true contact area is illustrated. To do so, the contact problems defined by cases of study I and II are solved under several configurations of the approach, in which φmaxhas been varied, keeping Luni<Lmax(setting 3 in Section 3.1). Figures 10e, f and 12b show examples of the resulting contact area and pressure element mesh that have been obtained for cases of study I and II under this setting of the approach.

The contact pressure distributions along the principal axes of the contact area of the solutions shown in Figure 10df are shown in Figure 11b. The contact pressure distributions along the principal axes of the contact area of the solutions shown in Figure 12a, b are shown in Figure 13. It both cases, it can be observed that increasing the value selected for φmaximplies 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 φmax, since the same values are obtained regardless of the value selected for this parameter. This is because when φmax<1, the accuracy in which the border of the contact area is computed depends only on the value selected for Lmax, as stated in Section 3.2.

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

5. Conclusions

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

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

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

  2. When Luni<Lmaxand φmax=0, adaptive mesh refinement is performed outside the true contact area. Under these circumstances, it can be observed that the computational cost of the approach is reduced by maximizing LmaxLini, while the accuracy of the solution remains unaffected.

  3. In last place, when Luni<Lmaxand φ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 φmaxis 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]).

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Victor Roda-Casanova and Francisco Sanchez-Marin (December 20th 2017). Adaptive Mesh Refinement Strategy for the Semi-Analytical Solution of Frictionless Elastic Contact Problems, Contact and Fracture Mechanics, Pranav H. Darji and Veera P. Darji, IntechOpen, DOI: 10.5772/intechopen.72422. Available from:

chapter statistics

214total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Improving Contact Load-Bearing Resistance of Ultrafine- Grained Materials Through Multilayering and Grading

By Dengke Yang, Jiangting Wang, Huimin Yang and Peter Hodgson

Related Book

First chapter

Manipulation of Tribological Properties of Metals by Ultrashort Pulsed Laser Micro-/Nanostructuring

By Quan-Zhong Zhao and Zhuo Wang

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us