Extending the BEM for Elastic Contact Problems Beyond the Half-Space Approach

The boundary element method (BEM) is widely used in fast numerical solvers for concentrated elastic contact problems arising from the wheel-rail contact in the railway industry. In this paper we extend the range of applicability of BEM by computing the influence coefficients (ICs) numerically. These ICs represent the Green’s function of the problem, i.e. the surface deformation due to unit loads. They are not analytically available when the half-space is invalid, for instance in conformal contact. An elastic model is proposed to compute these ICs numerically, by the finite element method (FEM). We present a detailed investigation to find proper strategies of FEM meshing and element types, considering accuracy and computational cost. Moreover, the effects of computed ICs to contact solutions are examined for a Cattaneo shift contact problem. The work in this paper provides a guidance to study fast solvers for the conformal contact.


Introduction
Contact problems have a broad range of applicability in engineering fields, such as wheel-rail contact in railway industry [15], investigation of friction and wear [5,26], rolling fatigue [7], ball bearing [9], etc. Such problems consider elastic contact between two bodies. When they are pressed together, an elastic field of stress, strain and displacement arises in each body. A contact area is formed where the two body surfaces coincide with each other. The stresses exerted on each body consist of normal stress and tangential stress. In the socalled normal contact problems, we solve for the contact area and the normal stress on it. When the two bodies are brought into relative motion and the relative velocity of the two surfaces is small, a creeping motion could be observed which is largely caused by the elastic deformation at the contact patch. In those parts of the contact area where the tangential stress is small, the surfaces of the two bodies stick to each other. Otherwise, local relative sliding may occur. This is studied in so called tangential contact problems, where the research question is to find the adhesion and slip areas, and the distribution of tangential stresses.
Solving contact problems has been a continuous research topic, where Johnson [12] and Kalker [13] contributed to three-dimensional contact problems. For more introduction of modern development, we refer to [25]. Numerical approaches for contact problems have been developed, typically in two classifications. One is the finite element method (FEM) [14,17,38], which is widely used, especially for problems involving large deformation and nonlinear materials. The other is the boundary element method (BEM) [1,13,22], which is applied for homogeneous material, under the assumption of small deformation and deformation gradients. This method particularly suits the concentrated contact problem well.

Strategies for solving concentrated contact problems
In concentrated contact problems, the contact area is considered flat since it is small compared to the dimensions of the contacting bodies. Examples can be found in Figures 1(a) and (b), where we depict the contact between the wheel tread and rail crown, and between a wheel flange and rail gauge face. The half-space approach is used, which approximates the contacting bodies as two semi-infinite solids bounded by the contacting plane [13]. The BEM transforms the three-dimensional (3D) boundary value problem to a two-dimensional (2D) boundary integral problem. The integral gives the relation between displacements and tractions. In the half-space approach, the corresponding Green's function is based on a solution by Boussinesq and Cerruti (see [12]). It is the so-called influence function, which expresses the displacements at one point caused by a load at another point on the half-space.
Solving such problems starts with the definition of the potential contact area, which contains the true contact region. Usually rectangular elements are placed in this 2D region, and piecewise constant functions are used to represent the tractions to be solved. This results in influence coefficients (ICs) from the influence function. Analytic coefficients are available [13], and have been incorporated with solvers for concentrated contact, e.g. [13,33,41].

Strategies for solving conformal contact problems
Conformal contact occurs when the rail and wheel are seriously worn so that a curved and much larger surface comes into contact, as seen in Figure 1(c). Figure 1. Two types of contact problems in wheel-rail contact [19]: concentrated contact in (a) and (b) with "flat contact areas", and conformal contact in (c) with a "curved contact area".
Other applications can be found, for instance, journal-bearing contact [36], microcontact printing [23], etc. The half-space approach is prohibited, since its assumptions of a flat and small contact area are not valid in these cases. The FEM is usually employed for such problems, e.g. in [6,20,31], due to its general applicability. However, since discretization by FEM covers whole contacting bodies, this method can be computationally slow and memory intensive. Moreover, focusing on the contact region, FEM often employs coarser meshes than BEM, which implies less detailed information in the region of interest.
One idea of fast solvers for conformal contact problems here is to apply the BEM, which means to only discretise the contact surface, and work on the surface integral similar as in concentrated contact [19,35]. Those fast solvers for concentrated contact can then also be applied. A difficulty is that the influence function (and further the influence coefficients) is not theoretically known in this case.
Therefore, the idea is to compute the influence function numerically by FEM, resulting in ICs. This approach has been applied in lubrication problems [8,16,37]. However, in the train's wheel-rail contact, i.e. a dry contact, different strategies of FEM meshing and element types are required. Furthermore, the effect of these computed ICs to final solutions is also different. Studying these two issues is an essential step for developing fast solvers for conformal contact in the railway industry. This approach was developed by Li in [19], where the rail was approximated by a quasi-quarter space with a distributed load and the FEM was applied. However, that work did not present detailed strategies for the mesh or element types. Moreover, there was no discussion of the error propagation from ICs to contact solutions. The approach was extended by Vollebregt [35], who computed influence coefficients on a 3D block caused by a piecewise constant load. He presented an initial meshing strategy of hexahedron elements, and a corresponding error propagation. The research is continued in the present paper. We will provide a detailed guidance to determine the mesh parameters for the hexahedron elements, and additional results on the error propagation. We will also include new results by adding a bilinear load. Another contribution is that prismatic elements will be investigated.

Content and structure of this paper
This paper serves as a preparation for the research of fast numerical solvers for conformal wheel-rail contact problems in the railway industry. In our earlier work, we have developed efficient solvers for concentrated contact problems, see [33,40,41]. Inspired by the use of the fast Fourier transform (FFT) technique in [21,24,28], our solvers employ FFT techniques, and show excellent efficiency from the point of computational time. This motivated us to extend the range of applicability of these solvers to conformal contact problems.
We investigate the influence coefficients in a concentrated contact setting, where analytic influence coefficients are known when a piecewise constant approximation is used as discretization of the tractions [12,13]. Secondly, there are analytic solutions such as the Hertz theory for normal contact [10] and the Cattaneo theory for the tangential shift problem [3]. The corresponding BEM results can provide a reference for studying conformal contact problems, for which there are no analytic influence coefficients, neither analytic contact solutions.
For conformal contact seen in Figure 1(c), one solves the displacements on the curved contact region of the rail. The FEM is then employed. Prismatic elements are generally applicable with their triangular faces on the curved boundaries of the cross section of the rail. Regarding such application, prismatic elements are also investigated for the model in this paper, which solves the displacements on a half-space.
In the model a unit load is applied on the region of one element, which brings a jump at the edge of this element to its neighbors. This implies a significant discontinuity in the traction boundary condition for the model as well, and may disturb the FEM solution. An alternative is to impose a bilinear load, which equals a "tent" shape at this region. Analytic investigation for this test setting can be found in [4], however, the implicit expressions for influence coefficients are difficult to process directly. A closed form solution is given in [18], however, on a triangular mesh. This is not useful for our work, which mainly focuses on rectangular meshes since these result in matrices of Toeplitz 1 structure so that the fast Fourier transform (FFT) can be applied to accelerate the corresponding matrix-vector products [33]. Therefore, numerical influence coefficients computed with a bilinear load on rectangular meshes are of interest.
The structure of this paper is as follows. In Section 2, an elastic model is built to compute influence coefficients. The results of this model are discussed in Section 3 with three different FEM mesh settings, from which insight is provided on the choice of meshing and element types, in view of accuracy and computational cost. In Section 4, a specific contact problem is solved, with numerical influence coefficients obtained by the recommended strategies. Error propagation is investigated separately in normal and tangential contact problems. The last section concludes this paper.

Influence coefficients on half-space
Influence coefficients are equivalent to displacements caused by a specific load. Evaluating the displacements under a given load is a basic problem in elasticity. It can be done by solving three elastic equations with some boundary conditions. These equations, in a quasi-static case, read [11]: 2. Strain-displacement equation: 3. Stress-strain equation (constitutive equation, or Hooke's law): In the above equations, subscripts i, j, k, l are related to the spatial coordinates x, y, z. The Einstein summation convention is used when repeated subscripts occur. Moreover, a comma means differentiation (i.e. σ ij,j = ∂σij ∂j ). Symbols σ, ε, F, u denote stress, strain, body force per volume, and displacement, respectively.
In the case of concentrated contact problems, an alternative way can be used to obtain displacements. It is based on four assumptions. First of all, the contacting region should be small compared to the dimensions of the contacting bodies. Secondly, these contacting bodies should be of homogenous linear elastic material. The third is small displacements and displacement gradients. Fourthly, inertial effects are ignored, which is a quasi-static case.
These assumptions allow for the use of the half-space approach, which approximates the two contacting bodies by two semi-infinite elastic bodies, bounded by the contacting plane. Usually a Cartesian coordinate system (O; x, y, z) is placed, with the origin at the center of the contact region, and the z-direction pointing to the outward normal of one contacting body. A halfspace is at one side of the bounding plane and can be denoted by, for instance, We only consider one half-space here. Displacements at one point on this half-space caused by a point load at another point were given as a function of relative distance between these two points by Boussinesq and Cerruti (see [12]). For instance, normal displacement u z at surface point (x, y) resulting from a normal point load P at another surface point (ξ, η) is given by [13]: where G is shear modulus and ν is Poisson's ratio. ρ = (x − ξ) 2 + (y − η) 2 is the distance between points (x, y) and (ξ, η). We drop the z-component for point coordinates here since all points of interest are on the surface of the half-space, i.e. z = 0. This equation has a singularity point when ρ = 0, at the point where the load is applied. This is not a serious problem since in practice the load is often distributed on a finite region [2]. The displacement caused by a distributed load can be obtained by integration of the displacement by the point load. Replacing normal load P by p z (ξ, η)dξdη in equation (2.4) and integrating on a certain contact area C, we obtain the relation between normal displacement and pressure (normal traction): where the kernel is called the influence function. It describes the influence of a unit normal load at point (ξ, η) on the normal displacement at point (x, y). The value of A zz (x, y, ξ, η) depends on the relative distance of these two points, rather than their own locations. This means that: Discretization is done by placing n := n x ×n y rectangular elements on a potential contact area which contains the real contact region. The corresponding element sizes are δx and δy, respectively. A cell-centered mesh arrangement is chosen. Using I, J ≤ n as element indices, a widely used approximation for pressure p z is by piecewise constant function f J (ξ, η): where f J (ξ, η) = 1 for (ξ, η) on element J and f J (ξ, η) = 0 elsewhere. Plugging equation (2.7) into integral (2.5), we arrive at: where e J denotes element J. Eq. (2.8) leads to a linear relationship: where matrix A zz contains all the influence coefficients (ICs) A zz IJ . Each of them is defined by: which results from integrating equation (2.5) on a single element e J w.r.t. an observation on element I. It indicates the physical meaning of the IC: A zz IJ is the normal displacement on element I, caused by a load f J on element J. In this case, the applied (normal) load f J is a unit piecewise constant function, and analytic solutions of integral (2.9) can be found in [13]. 2 Moreover, Eq. (2.9) also implies that the value of the IC is related to the area of element e J and material parameters G, ν.
A similar derivation can be done for other ICs, such as A yx , the y-component of displacement caused by a unit load in the x-direction.

Model
The model for computing the ICs is based on the contact problem to be solved. Assume that the contact problem considers two bodies of homogeneous linearly elastic material, with shear modulus G and Poisson's ratio ν. Recall that n = n x × n y rectangular elements are placed in the potential contact area, with corresponding mesh sizes δx and δy.
This model is inspired by the fact that an IC is equivalent to the displacement at one surface point resulting from a unit load at another surface point. We have to deal with two issues. One is to approximate the semi-infinite halfspace, and the other is to provide a proper boundary value problem. To approximate the half-space which represents one contacting body, we consider a block of size Figure 2(a)). This block is of the same material as the contacting body, and L → ∞ gives the half-space model. A unit load is applied at the center of the top surface, and we solve for the resulting displacements on this surface, which are the corresponding ICs. Due to the symmetry and anti-symmetry of the displacements, the computational domain can be reduced to a quarter of the block, [−L, 0] × [−L, 0] × [−L, 0], which is shown in Figure 2(b). Table 1. Boundary conditions at the front surface (x = 0), at the right surface (y = 0) and at the top surface (z = 0) for different ICs. On the other three surfaces the boundary condition is given by u = 0.

ICs
Front surface (x = 0) Governing equations for this model are the basic equations in elasticity, as given by equations (2.1)-(2.3). Boundary conditions can be different for different ICs. Define displacement by u = (u x , u y , u z ) T , and stress by σ = [σ x , σ y , σ z ] T , then Table 1 gives boundary conditions at the front, right and top surfaces of the block. At the other three surfaces the boundary condition is given by u = 0.
We explain the boundary conditions on the top surface in Table 1. The top surface of the block is shown in Figure 3. A 2 × 2 grid with h c = δx = δy is used for the contact problem. Due to the cell-centered mesh, the ICs for the center points of these four elements need to be computed. Since the value of the IC depends only on the relative distance of a pair of points, we can place the right-bottom element of the contact grid centered at the origin and compute displacements for these four points w.r.t. the origin. Moreover, a piecewise constant representation of tractions is employed, which implies that a unit load is applied on the region of the center element (the shaded area in Figure 3). Since we only consider a quarter of the block, the load in the model is only given on a quarter of the center element.
With governing equations (2.1)-(2.3) and boundary conditions, we can employ the finite element method (FEM) to solve for the displacements. In Figure 3, differently from the 2 × 2 contact grid, other meshes can be used when solving this model by the FEM. (Be aware that the meshes for the contact problem, and for the FEM model are different.) A linear interpolation technique is applied to convert the resulting displacements to the center points of the contact mesh, which leads to the required coefficients. Remark 1. When using a bilinear approximation of pressures, i.e. the applied load is bilinear in the model, boundary conditions on the top surface z = 0 are different. For instance, for the IC caused by a normal load, i.e. A xz , A yz , A zz , the condition reads [32]: For the IC caused by a load in x-or y-direction, one needs to replace the left-hand side of this condition by σ x (x, y, 0) or σ y (x, y, 0), respectively.

Influence coefficients A zz by a piecewise constant load
We know the analytic ICs when pressures are represented by piecewise constant functions [12,13]. In this section, the model proposed in Section 2 is used to compute the IC A zz , which is the normal displacement caused by a unit normal load. The resulting IC value can be compared with the analytic IC. Tests in this section aim at providing insight on mesh strategies and element types, considering the accuracy and computational cost. The contact problem which requires the numerical IC is defined as follows. Consider a sphere pressed on a plane. Both bodies are of the same material, with shear modulus and Poisson's ratio G = 82000 N/mm 2 , ν = 0.28, respectively. In this case, the analytic IC divided by 2 is the displacement on one of the bodies. The potential contact area is [−4, 4] × [−4, 4] mm 2 . On this area, four grids are placed: 10 × 10, 20 × 20, 40 × 40, 80 × 80, with mesh sizes h c = 0.8, 0.4, 0.2, 0.1, respectively.
As a reminder, subscript c refers to the contact problem, subscript s in the following refers to the model, which is solved by finite element package SEPRAN [27]. For the solution of the FEM discretization, a conjugate gradient method with an ILU preconditioner is employed with stopping criterion where u is displacement and the superscripts are iteration indices. To compare with the analytic IC on a certain contact grid, linear interpolation is used. We define our target of approximation by: where A zz 0 is the analytic coefficient (divided by 2) at the origin, and its approximation by the FEM is F zz 0 . We examine the errors at the origin since here the largest error occurs. The meshes that are recommended based on target (3.1) result in satisfactory contact solutions (see Section 4).

Mesh strategies
Remember that our FEM model needs to approximate the semi-infinite halfspace. Hence, the FEM mesh should be easily extended to a very large block without involving too many elements. Since the displacements decrease outside the loaded region, more elements can be used near the origin, and much fewer elements far away from it. Regarding the type of elements, we first use the hexahedron element which is easy to be generated for a 3D block. The corresponding mesh strategy is called "Mesh1", shown in Figure 4 Outside this cube, the element size increases by a factor f leaving the loaded area. Moreover, edges that do not contain the small cube are divided uniformly into m intervals. In this mesh, stretching factor f , cube size l and its element size h s are to be determined, so that accurate numerical ICs result.
Another employed element is the prismatic element. As mentioned in Section 1, in conformal contact, one computes ICs considering a load on the curved contact surface of the rail, see Figure 1(c). The cross section of the rail has curved boundaries, and it is well fitted by triangles. Moreover, the rail can be regarded as infinitely long, so rectangles are used for the other surfaces. This yields the prismatic elements.
Assuming that the front and back surfaces of the 3D block correspond to the cross section of the rail in Figure 1(c), these two surfaces should consist of the triangular faces of the prismatic elements, as shown in Figure 4 This prismatic mesh is called "Mesh2", which also has the small cube configuration with a uniform mesh to cover the loaded region. The mesh strategy for the front and back surfaces is the same as the front surface of hexahedron Mesh1, except for using triangles instead of rectangles. The top surface can be divided into four parts. One corresponds to the small cube. The second is a rectangle along the x-direction, with its width equal to cube size l. Its element size, along the x-direction, increases by a factor f outside the cube, but the element size in the y-direction is always h s . This will result in extremely stretched elements when the block size L gets large. The rectangle along the y-direction has similar structure. In the fourth part of the top surface, the element size increases in both directions. The right surface has a similar structure.
Next, we will test different mesh parameters on Mesh1 and Mesh2. Three tests are defined in Table 2, w.r.t. h c , the element size of the grid for the contact problem. TESTs 2 and 3 are based on a smaller cube near the origin than TEST 1. Moreover, TEST 3 has a smaller stretching factor f than TEST 2.

A zz on hexahedron Mesh1
We start by approximating the IC A zz for the 10 × 10 contact grid. Figure 5 shows results for TEST 1 using linear and quadratic elements, as block size L increases. Linear elements fail to achieve the target accuracy (3.1), denoted by the dashed line. However, quadratic elements exhibit very satisfactory convergence. Therefore, quadratic elements will be employed in the subsequent tests. Focusing on quadratic elements, the corresponding errors decrease as L is doubled, and they stabilize for L ≥ 512. The reason for this is that the error consists of two parts. One is the error of the continuous model, where an infinite domain is approximated using a finite block size L. The other part is the finite element discretization, approximating the continuous solution on the finite domain. When L is small, both types of error occur and the first part is dominating. When L is larger, e.g. L ≥ 512, the first part is annihilated and only the second part remains. Then it doesn't help anymore to increase L further. TEST 2 and TEST 3 use a smaller uniform cube, and the latter has a smaller stretching factor than the former. Figure 6 shows the results of these two tests. Since TEST 3 employs a finer mesh, it results in smaller errors than TEST 2. At the same time, it requires more elements and thus CPU time. But the cost is acceptable. Comparing with TEST 1, we observe that TEST 3 can achieve similar errors, however at a lower cost (since it has a coarser mesh).
Moreover, the CPU time of TEST 3 shows an almost linear dependence on L in Figure 6(c). This is due to the strategy to generate hexahedron Mesh1, where the number of elements n is approximately linearly dependent on L, as can be seen in Figure 6(b). The FEM discretization results in a linear system with a sparse coefficient matrix, where the number of non-zero entries in each row is much smaller than n. Remember that a conjugate gradient (CG) method with an ILU preconditioner is applied to solve this system. The results in Figure 6 (and Figure 7) indicate that the complexity of this algorithm is O(n 7/6 ) (consistent with the statement of [29]). For such a sparse matrix, the total number of CG iterations is mildly depending on L. For instance, there are 76, 78, 82, 85 CG iterations when L = 128, 256, 512, 1024, respectively.
We need to balance the accuracy and the computational cost. On the one hand, a coarser grid yields cheaper solutions, for instance, using larger stretching factors such as f = 2 and f = 2.5, or a smaller uniform cube l/h c = 0.5. However, the corresponding accuracy may be unsatisfactory. A finer grid gives rise to a better accuracy but is more expensive. Figure 7 gives the results using TEST 3 on the 20 × 20, 40 × 40 and 80 × 80 grids with h c halved. It can be seen from Figure 7 reach the target accuracy when L/h c ≥ 80 is satisfied. Moreover, when L/h c is the same for these three grids, the corresponding errors are also similar. In other words, by shifting the line with red circles to the right, it will coincide with the line with black stars. The same observation can be found in Figures 7(b) and (c). The reason is that the parameters are determined based on the scaling factors w.r.t. the mesh size of the contact problem h c . When L is fixed, a finer contact grid requires computing using a finer FEM mesh. Therefore, the corresponding error is smaller than that on a coarser contact grid, as seen in Figure 7(a). We only show results with prismatic elements for a 10 × 10 contact problem in Figure 8. It can be found that TEST 3 leads to the smallest errors when L increases. However, the corresponding cost is the highest (the CPU time with L = 1024 is so high that we do not show it in this figure). The errors by TEST 1 and TEST 2 are comparable. Since the latter requires fewer elements and less CPU time, it is the better choice for prismatic elements on Mesh2.

A zz on prismatic Mesh2
A difference with the behavior for the hexahedron mesh observed in Figure 6 is that the errors at the origin with prismatic elements are decreasing less regularly when L increases, and stagnate at a higher level than when hexahedron elements are used. This is attributed to the large stretching of the prismatic elements when a large L value is used, which yields an ill-conditioned matrix from the FEM discretization, and hence the results are spoiled to some extent.

Discussion on mesh strategies
Based on the above results, we conclude regarding the choices of the mesh parameters in Table 3 for hexahedron Mesh1 and prismatic Mesh2. With these parameters, target accuracy (3.1) can be achieved. We can see that the same ranges apply for domain size L, cube size l and cube element size h s between the hexahedron and the prismatic meshes. Stretching factor f is optimal in different ranges. The optimal choices for Mesh1 and Mesh2 are TEST 3 and TEST 2, respectively (See Table 2 for the definition of these two tests). Other conclusions are as follows: • Quadratic elements are preferable on both Mesh1 and Mesh2.
• The errors become stable when L/h c ≥ 640 on Mesh1, and when L/h c ≥ 160 on Mesh2.
• The same mesh strategies can be used to compute other ICs such as A xx , A yy , and A xy , using the corresponding boundary conditions as defined in Table 1. Similar results are found for the numerical accuracy and the computational time.
• The same mesh strategies are also applicable to compute ICs when prescribing a unit bilinear load.
• Smaller errors can be attained with Mesh1 than with Mesh2, as described in Section 3.3. A proper value of domain size L needs to be chosen for the prismatic Mesh2.

Error propagation in a Cattaneo shift problem
Errors in influence coefficients (ICs) can disturb the final solution of a contact problem. To investigate this, a Cattaneo shift contact problem [3] is solved, which involves a 3D frictional contact. It considers a sphere which is first pressed and then tangentially shifted on a plane, so that partial slip on the contact area may occur. The test problem was specified in Section 5.2.1.1 of [13] and in Section 5.1 of [34]. Both sphere and plane are of the same elastic material, a quasi-identical case. In this case, the problem can be decoupled into a normal and tangential contact problem, since the results of the latter do not influence the former [13]. We solve these two subproblems separately in this section. Analytic and numerical ICs are used in the solution procedure, and the contact solutions are compared.

The normal contact problem
The normal contact problem is based on a normal force applied to the sphere on the plane so that a contact area is formed due to the deformation. The sphere and plane have the same shear modulus and Poisson's ratio, G = 200 N/mm 2 and ν = 0.42, respectively. The radii of the sphere are R x = R y = 50 mm. The contact area is a circle with radius 1 mm. According to the Hertz theory [12], the normal force is F z = 9.1954 N , the approach, i.e., the maximal penetration, is δ = 0.02 mm, and the maximal pressure is p max = 4.3905 N .
In the numerical implementation, approach δ = 0.02 mm is prescribed. The contact area and pressure on it are solved for. We show the results for the normal force F z , which is the integral of the pressure on the contact area. Moreover, the errors in the maximal pressure p max and in the contact area are also analyzed. The potential contact area is set to [−1.2, 1.2]×[−1.2, 1.2] mm 2 . The undeformed distance function reads h(x, y) = 1 2Rx x 2 + 1 2Ry y 2 − δ. The algorithm NORM [13] is used, with stopping criterion chosen as where d is the defect of governing equation A zz p z + h = 0, and the root-mean square norm is used.
Only IC A zz is required in this normal contact problem. Both analytic and numerical ICs based on our model are employed in the numerical experiment. Relative errors w.r.t. the Hertzian solutions are shown. For example, the error in the normal force is defined by |Fz−Fz| Fz , whereF z and F z are forces provided by the numerical method and analytic solution, respectively.
We firstly examine the IC obtained by a unit piecewise constant load, then the ICs by a bilinear load. The results are also compared.

Analytic ICs
First of all, we use the analytic ICs in this normal contact problem. Relative errors in the normal force F z , in the maximal pressure p max and in the area of contact region are presented in Figure 9. It is found that the errors in the normal force are already quite small. Errors in the maximal pressure show a linear reduction when the grid is refined. The errors in the contact area are the largest, which can be explained by the addition or removal of complete elements in the numerical treatment where the contact area changes.

ICs by a piecewise constant load on hexahedron Mesh1
In this experiment numerically computed ICs on Mesh1 are used. Figure 10(a) shows relative errors in the resulting normal force for three different discretizations of the contact problems: 10 × 10, 40 × 40 and 160 × 160 grids. It can be seen that these errors decrease as domain size L increases on the three grids. Moreover, the errors using 40 × 40 and 160 × 160 grids are very similar, and they are smaller than using the 10 × 10 grid when L increases.   Figure 10(b) presents the relative errors in the maximal pressure p max . As L increases, the error for the 10 × 10 grid decreases until it stagnates, then it increases and stabilizes. Such stagnation also occurs with the errors related to the 40 × 40 and 160 × 160 grids (the latter is not clearly shown since stagnation is at L = 1024. If L is larger, the curve will increase again). However, it is difficult to know beforehand the value of L for which the maximal pressure stagnates. Moreover, when L reaches the largest value, i.e. L = 1024, the finer the grid is, the smaller the resulting error.
Errors in the contact area are presented in Figure 10(c), where the finest grid always yields the smallest errors for different values of L. Moreover, as L increases, errors on these three grids reduce and then stagnate.
Based on the above discussion, we come to the conclusion that it is optimal to choose L as large as possible on this mesh.

ICs by a piecewise constant load on prismatic Mesh2
ICs computed on prismatic Mesh2 are employed here. In Figure 11(a), errors decrease as L increases, and then they increase when L gets larger. Using 40 × 40 and 160 × 160 grids results in smaller errors than the 10 × 10 grid. Moreover, these two fine grids exhibit the same errors. This situation is the same as when using the IC on Mesh1 in Figure 10(a).   Figure 11(b) presents a stagnation in each curve for the error in the maximal pressure p max . Errors in the contact area in Figure 11(c) indicate that for a fixed value of L, the finer the discretization, the smaller the error will be.
Notice that in Figure 11, the results on 40×40 and 160×160 grids are lacking when L gets very large. The reason is that the ICs are not computed in this case. In fact, highly distorted elements occur on Mesh2. As mentioned in Section 3.3, large L values give rise to over-stretched elements and ill-conditioned FEM matrices. Therefore, it is important to carefully find a proper value for L on Mesh2, neither too large nor too small.

ICs by a bilinear load
Here, the ICs have been computed also with the bilinear load as described in Remark 1. Since the analytic values cannot be accessed easily, no comparison was made in Section 3 between analytic ICs and numerical ICs. However, it is possible to use these ICs in the contact model and study the propagation of errors, as the above discussion for the ICs computed using the piecewise constant load. The resulting figures are very similar to Figures 10 and 11, and are therefore omitted for brevity.

Discussion on the normal contact problem
We compare the influence of all ICs on the normal contact problem. Numerical ICs are obtained on Mesh1 with L = 1024 and on Mesh2 with L = 128. Relative errors in the normal force are shown in Figure 12(a). We find that, on the one hand, both piecewise constant and bilinear loads result in similar errors on a very fine discretization. On the other hand, errors on Mesh1 are smaller than those on Mesh2. For instance, on the finest 160 × 160 grid, the former reaches a relative error of 0.1% and the latter around 0.8%. Moreover, these numerical ICs cannot yield normal forces as accurate as the analytic ICs, which results in errors smaller than 10 −6 .   Figure 12(c), where we notice that the ICs based on a bilinear load on Mesh2 can even result in better accuracy than the analytic ICs.

The tangential contact problem
In the tangential contact problem, the contacting sphere is shifted in the xdirection by a tangential force F x = (7/8)µF z . Here, µ = 0.4 is the friction coefficient of Coulomb's frictional law. The sphere sticks to the plane where the tangential traction is small. Local sliding occurs when the traction reaches the frictional traction bound. According to Cattaneo's theory [3,13,34], the adhesion area is a circle with radius θ = 0.5. In this case, the rigid shift of the sphere is given by w x = 0.00817 mm.
In the implementation, we have prescribed the rigid shift w x = 0.00817 mm. The tangential traction, the adhesion and slip areas are to be solved. Integra-tion of the tangential traction over the contact area yields tangential forces F x , F y . The required ICs are A xx , A yx , A yy , A xy , where A yx = (A xy ) T . The solver is the TangCG algorithm, proposed in [41]. The iterations of TangCG terminate when ||p k+1 t −p k t || ||p k+1 t || < 10 −6 , where p t is the tangential traction consisting of x-and y-components. The solution procedure starts with the computed contact area and the corresponding pressures resulting from the normal contact problem where the same strategies are used to compute the ICs for both the normal and tangential contact problems. This means that, for instance, if A zz for the normal contact problem is obtained by a bilinear load on Mesh1, this strategy is also employed to compute A xx , A xy , A yy for the tangential contact problem.
We only compare tangential force F x in this case. The analytic Cattaneo solution for F x is not completely relevant here, since it does not involve the force in the lateral direction w.r.t. the shift direction. Therefore, rather than using this solution as the reference for comparison, we use F x = 3.2216 N , which is obtained by the TangCG algorithm with analytic ICs when h c → 0.  First, we use the ICs computed by a unit constant load on hexahedron Mesh1, and the resulting errors in tangential force F x are shown in Figure 13(a). These errors decrease as domain size L increases. Errors on the 10 × 10 grid sometimes are smaller than those on finer grids. This may due to the fact that a finer grid requires a larger number of ICs and hence errors may accumulate in the evaluation. Figure 13(b) presents errors when the ICs are computed on prismatic Mesh2. As the domain size increases, these errors first are reduced, and then increase again. The same observation can be found when the ICs are obtained by a bilinear load, so we do not show the corresponding results here.
To compare these ICs, again we choose L = 1024 for Mesh1 and L = 128 for Mesh2. A piecewise constant load as well as a bilinear load are used. Figure 14 shows the errors in tangential force F x with the resulting ICs. Using analytic IC yields a linear reduction of errors, and the numerical ICs obtained by hexahedron Mesh1 give rise to smaller errors than those by prismatic Mesh2. For instance, the former achieves around 0.03% and the latter possesses 0.17% on the 160 × 160 grid. The piecewise constant load and bilinear load yield similar results. However, we conclude that "engineering accuracy" can also be Here, "PC" and "Bi" denote piecewise constant and bilinear loads, respectively. The horizontal axis shows the number of elements in each direction of contact grid nc × nc.
achieved on the prismatic Mesh2.

Conclusions
In this paper, we studied the numerical calculation of the influence coefficients (ICs) that are used in fast solvers for concentrated contact problems. These ICs are displacements caused by a unit load, either a piecewise constant load or a bilinear load, depending on the representation of the unknown tractions. To compute these ICs, an elastic model was provided, including governing equations and boundary conditions for different ICs. The finite element method is employed via package SEPRAN. Considering the accuracy and computational cost, we recommend hexahedron Mesh1 and prismatic Mesh2, with their own specific mesh parameters after several numerical tests. The former is easily extended to a large computational domain that approximates the half-space. The latter is more generally applicable, while a proper domain size needs to be chosen carefully in order to prevent over-stretched elements.
The numerically computed ICs are incorporated into the solution procedure of a Cattaneo shift problem, which is decoupled into a normal and a tangential contact problem. The effects of these ICs on the contact solutions are discussed. ICs obtained on hexahedron Mesh1 show better performance than those on prismatic Mesh2. However, the prismatic mesh is favorable for engineering applications. The performance does not differ much between ICs that are computed by a piecewise constant load or a bilinear load. This means that using piecewise constant and bilinear representations for unknowns result in the same accuracy, in agreement with the corresponding statements in [13,32].
A different strategy for computing ICs would have been to use the adaptive finite element method with a-posteriori error estimation (see e.g. [30,39,42]). This will refine the mesh automatically in the vicinity of the loaded region, where the most rapid variation of the solution occurs. The goal of such auto-matic refinement is to provide meshes that are quasi-optimal in some sense. The benefit of designing the meshing strategy by hand is that the precise demands of the application can be accounted for, such as the regular grid structure at which the ICs are required. An interesting topic for further study is therefore to compare the performance of adaptive FEM to the methods that are presented here.
The work in this paper provides a guidance for developing fast solvers for conformal contact problems. Similar model and meshing strategies can be used there to compute ICs on curved domains. Similar errors in the ICs and similar error propagation may be expected, such that meshing parameters may be chosen along the lines of the results presented here.