Multigrid numerical simulation of contact mechanics of elastic materials with 3D heterogeneous subsurface topology

Tribology International 92 (2015) 233–245

Contents lists available at ScienceDirect

Tribology International journal homepage: www.elsevier.com/locate/triboint

Multigrid numerical simulation of contact mechanics of elastic materials with 3D heterogeneous subsurface topology Hugo Boffy n, Cornelis H. Venner Faculty of Engineering Technology, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands

art ic l e i nf o

a b s t r a c t

Article history: Received 29 January 2015 Received in revised form 2 June 2015 Accepted 15 June 2015 Available online 23 June 2015

Contact phenomena between deformable bodies are a common problem in engineering. The surface stress distribution and subsurface stresses depend on many parameters. In view of increasingly strict tolerances the effects of local variations in material properties need to be accurately predicted. In this paper a multigrid solution method is presented for contact problems between three dimensional elastic heterogeneous materials. The contact problem is incorporated as boundary condition in the multigrid solution of the displacement equations for the volume. First, validation results are presented. Subsequently a study is presented for soft and hard clusters of inclusions. Finally, results are presented for a contact problem involving a realistic case of a polycrystalline material representative for applications with ceramic materials. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Multigrid Linear elasticity Heterogeneous properties Contact mechanics

1. Introduction Prediction of the behavior of dry and/or lubricated contacts is a crucial part of design of many forming processes, and in machines, and machine elements. Considering elements with contraform contacts such as in gears and bearings, failure is often due to subsurface fatigue which initiates in regions of concentrated stresses below the surface, induced by the contact normal stress (pressure) and shear stress distribution moving over the surface under rolling/sliding. The contact stress distribution depends on the geometries of the mating bodies, the microgeometry (surface topography, roughness), the loading and operating conditions, the material properties, including the variation of the material properties due to the structural variations of the material: grains, coatings, fibers. As a result of increasing demands on efficiency, reduction of material use, and the use of more complex and composite materials prediction of many small scale effects is becoming increasingly important. Contact problems are studied in many disciplines in engineering, for example using the Finite Element Method (FEM). Hughes et al. [1,2] worked on the formulation for compatible meshes in the case of contact-impact problems accounting for phenomena such as stick, slip or frictional sliding conditions. Curnier [3] constructed a theory of friction for small slip inspired by the theory of plasticity at small strains. Simo et al. [4,5] developed a segment approach to discretize the contact interface without

n

Corresponding author. Tel.: þ 31 53 489 24 88. E-mail address: [email protected] (H. Boffy).

http://dx.doi.org/10.1016/j.triboint.2015.06.015 0301-679X/& 2015 Elsevier Ltd. All rights reserved.

having the need for element nodes to match each other. Bathe and Chaudhary [6] also used the contact segments and developed a solution procedure for the analysis of planar and axisymmetric contact problems involving stick, frictional sliding and separation under large deformations. For more details about FEMs for contact problems, see, e.g., [7]. More recently, Schröder [8] developed a mixed FEM of higher-order to deal with idealized frictional problem. Hild and Renard [9] present an improved the FEM error analysis by using estimates on Poincaré constants. Coupled methods have also been developed to deal with elasto-plastic contact problems such as the Element-Free Galerkin-Finite Element (EFGFE) [10] or elastoplastic Finite Element and elastic Boundary Element Methods (FEM-BEM) [11]. The common approach to contact mechanics in tribology, where the contact region is often small relative to the radius of curvature of the contacting elements is to assume two bodies with a particular (local) undeformed geometry brought into contact and for the deformation to assume the bodies to be semi-infinite, homogeneous and elastic, so that the surface deformation can be expressed in terms of integral transforms. The contact pressure is then solved from the Fredholm integral equation of the first kind which equates the undeformed gap to the (normal) displacement, subject to the conditions of Hertz–Signorini–Moreau and the condition that the integral over the pressure field should balance the externally applied load. For paraboloid surfaces the solution is the well known analytical solution found by Hertz [12]. For general shapes a numerical approach is required for which Kalker et al. [13] have been some of the pioneers. For the numerical solution of the pressure a direct or an iterative method can be used. A complication arises from the

234

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

surface deformation integrals in the analytic formulation which appear as multisummations in the numerical model. They lead to a full matrix which is expensive to solve directly. Also when using an iterative solver, they have to be evaluated many times during the iterative process which also leads to excessive computing times. In the past decades various approaches have been presented for the faster evaluation of the integrals, and the solution of the problem. Brandt and Lubrecht introduced Multilevel Multi-Integration (MLMI) for the fast evaluation of the integral transforms integrated in a Multigrid solver with distributive relaxation [14] for the contact problem to obtain a solver with ðN log ðNÞÞ work count if N is the total number of unknowns, and used the same technique (MLMI) to obtain the subsurface stresses [15]. Polonsky and Keer [16] used the conjugate gradient method for the solution of the pressure field combined with Fast Fourier Techniques (FFT) for the fast evaluation of the deformation integrals. A comparison of FFT and Multigrid/MLMI can be found in [17]. The topic of solving the non-adhesive contact problem in this standard formulation remains to attract attention, e.g., see Vollebregt [18]. For the case of steady state contacts separated by a lubricant film the approach can be similar, except that the Reynolds equation is added to the system, e.g., see [19]. However, for transient cases, moving roughness, sliding surfaces the problem of elasto-hydrodynamic lubrication is quite different from the non-adhesive contact problem. Combined with today's generation of computers the aforementioned methods allow detailed simulation of contact problems in which, e.g., the effect of the surface micro-geometry is taken into account, including cases of multiple islands of contact. The computed pressure fields at the surface can be used as input to stress calculations to, e.g., determine the effects on contact fatigue in rolling bearing contacts [15]. It is well established that especially for rolling contacts at sufficiently high loads or in the starved regime, the predicted stress fields from a non-adhesive contact analysis are quite accurate. With load conditions getting more severe, and tolerances narrower, effects of variations in properties become more important, e.g., not only failure due to inclusions or (clusters) of pores needs to be predicted accurately, but also cases with materials which were so far successfully modeled as homogeneous, need to be reconsidered for the small scale topology variations due to granular structure. Also, failure criteria for specifically designed coated and composite materials, need to be designed by numerical simulations, of course exploiting observed scaling and self similarity so as to minimize the number of simulations that needs to be done and to derive models and function fits that correctly reflect the underlying physics. The interest in non-homogeneous materials is not new, modeling of layered surfaces and coatings was already done in the 1960s. An analytical solution for the case of a rigid spherical indenter contacting an elastic layer on a rigid substrate has been presented by Vorovich and Ustinov [20]. This approach was extended to an elastic indenter by Keer [21]. Numerical methods were developed to solve contact problems involving one or several layers in two and three dimensional problems: Kennedy and Ling [22] used the finite element method (FEM) to analyze the indentation of elastic and elastic-plastic material. Plumet and Dubourg [23] determined the contact parameters and the stress field for the case of a multilayered body, based on linear elasticity theory, integral transforms (IT), Fast Fourier Transform (FFT), and contact analysis with friction. Polonsky and Keer [24] used a CG/FFT method for layered contacts and to compensate for the periodicity error introduced by the FFT method, used a special correction term which is computed using the MLMI technique developed by Brandt and Lubrecht [14]. Cai and Bhushan [25] used an FFT based approach and accounted for surface roughness. Watremetz et al.

[26] and Boffy et al. [27] used a finite difference discretization of the equations of elasticity in the material volume and a multigrid solution method to solve for the displacement field and stresses in coatings and functionally graded materials subjected to a surface pressure distribution. This approach was used within the conjugate gradient method for the solution of the contact problem. Heterogeneity of the subsurface material may appear also in the form of defects or inclusions. Since the pioneering modeling work of Eshelby [28,29], the effects of single or multiple inclusions on the stress field have been studied by many authors. Recently, Nelias et al. [30,31] and Zhou et al. [32,33] developed numerical approaches to perform parametric studies of the influence of inclusions on the contact parameters. Nelias et al. use the Eshelby model and a 3D-FFT approach to solve the displacement and stress fields on the surface and in the subsurface. Zhou et al. [34] also used a semi-analytical model based on Eshelby's solution for the subsurface. They derived a (integral/multisummation) formulation for the surface deformation and solved the contact problem using the conjugate gradient method. A disadvantage of FFT based methods is that they tend to loose efficiency when the 3D surface topology becomes more and more irregular, e.g., arbitrarily distributed inclusions, grains structures, or irregular fibers, whereas such cases need to be considered, e.g., for ceramics and composite materials. Boffy and Venner [35] demonstrated that a multigrid algorithm for solving the 3D elastic problem does not have this disadvantage when properly developed. This approach can be used as an auxiliary “routine” inside any numerical contact pressure solver, e.g., within a conjugate gradient [36] solver. However, by their nature multigrid techniques allow a more efficient approach, namely simultaneous solution of the elastic problem with the contact problem incorporated in the boundary condition. In that case the entire problem is solved in an amount of work that is just a little larger than just solving the 3D elastic problem only once. This approach is offered in the present paper. In addition to the explanation of the methodology, results are shown for representative applications. First for simple validation cases of a homogeneous and a coated material. Next for material consisting of grains with interstitial glue characteristic for a ceramic application.

2. Theoretical model The problem considered is that of an indenter of given geometrical shape contacting a 3D volume of heterogeneous elastic material. In Section 2.1 the equations for the volume are given. In Section 2.2 the boundary conditions are described.

2.1. Equations The problem can be described as the solution of the unknown displacements u, v and w in a 3-D domain Ω from the Navier– Cauchy equations: ðλuj;j Þ;i þ ðμui;j Þ;j þ ðμuj;i Þ;j ¼ 0;

i; j ¼ 1; 2; 3:

ð1Þ

The indices 1, 2 and 3 refer to the x; y and z directions respectively with the displacements u1 ¼ u, u2 ¼ v, and u3 ¼ w. λ and μ are the Lamé coefficients which are assumed to vary as a function of space:

λðx; y; zÞ ¼

Eðx; y; zÞνðx; y; zÞ ; ð1 þ νðx; y; zÞÞð1  2νðx; y; zÞÞ

ð2Þ

μðx; y; zÞ ¼

Eðx; y; zÞ : 2ð1 þ νðx; y; zÞÞ

ð3Þ

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

Ω will be taken as a 3D cubic domain, which may be representative of a sub-volume of a heterogeneous material. 2.2. Boundary conditions Without loss of generality it is assumed that a rigid indenter deforms the surface z¼ 0 which is the top of the cubic domain. At the bottom surface the displacements are specified: u ¼ v ¼ w ¼ 0. The vertical faces of the cubic domain are assumed to be stress free. This choice in the boundary conditions is not a restriction, one can use for example the results (stress or displacement fields) obtained by a macro-scale simulation (using FEM, FFT or other methods). Here only the equations for the top surface are given. At this surface a zero shear stress will be assumed: ∂u ∂z

∂w ¼ 0; ∂x

ð4Þ

∂v ∂z

∂w ¼ 0: ∂y

ð5Þ

σ xz ¼ 0- þ σ yz ¼ 0- þ

Next the contact condition should apply. Let the gap between the indenter and the surface be given as hðx; yÞ ¼  δ þ gðx; yÞ  wðx; y; z ¼ 0Þ; with gðx; yÞ being the shape of the indenter, and approach. The contact equation is h¼0:  δ þ gðx; yÞ  wðx; y; z ¼ 0Þ ¼ 0:

ð6Þ

δ the mutual ð7Þ

The contact normal pressure is given by σ zz ¼  p leading to the condition:
 ∂w ∂u ∂v þ þ p ¼ 0; ð8Þ ðλ þ 2μÞ þ λ ∂z ∂x ∂y subject to p Z 0 (adhesion is neglected). If the boundary equations give a negative pressure this would imply that the gap is open, so Eq. (7) cancels, and Eq. (8) applies with p ¼0. So, the boundary conditions are either a set of four equations (4)–(8) for u; v; w and p, or a set of three (4), (5) and (8) for u; v and w with p¼ 0. The formal description of the normal contact boundary condition is the Hertz–Signorini–Moreau condition: p Z 0, h Z 0, and ph ¼ 0. Finally, the solution of the pressure at the top surface is subject to the condition of force balance: Z pðx; yÞ ¼ f ð9Þ z¼0

where f is the external load applied to the indenter. This equation determines the mutual approach δ in Eqs. (7) and (6).

3. Numerical solution The equations were discretized with second order accuracy on a uniform grid using a collocated approach. The motivation for this choice and the discrete equations are given in Appendix A. A Gauss–Seidel iterative algorithm to solve the equations is described in Appendix B. However, iterative schemes like this have a grid dependent (slow) convergence rate. Error components with a wavelength of the order of the domain are reduced only by 2 a factor 1  Oðh Þ if h is the mesh size. As a result solving the problem to an error below discretization error will cost Oððp=dÞNð2=d þ 1Þ log ðNÞÞ operations if N is the total number of gridpoints, p the order of discretization, and d the dimension of the problem. This will lead to excessive computing times for large N and seriously limits the amount of detail that can be described on a grid within acceptable CPU time which will significantly limit the capability to study the influence of variation of material properties on a granular scale. High frequency components with

235

a wavelength of the order of the mesh size on the other hand are very efficiently reduced, typically by an order of magnitude in only a few iterations. This behavior can be exploited to obtain a very efficient algorithm by using multigrid techniques, as introduced by Brandt [44]. The iterative process is used only a few times to reduce the high frequency error components and smooth the error. Subsequently a usually twice coarser grid is used to approximate and solve the remaining smooth error components. This result is used as a correction to the fine grid solution. This approach, applied recursively leads to a Multigrid coarse grid correction cycle which, properly designed, can have a grid independent convergence rate. In addition to using coarser grids for convergence acceleration they can also used to efficiently generate an accurate first approximation on the target grid. The resulting algorithm allows solving the problem up to the level of the discretization error in an amount of work O(N), i.e. proportional to the number of unknowns. As a result the method has the potential to use dense grids with many unknowns on small scale computers facilitating detailed description of a material topology. In this paper it is demonstrated that this prospect is realistic for 3D elastic contact problems. A flow diagram of a multigrid algorithm illustrating the use of the multiple grids is given in Fig. 1. An introduction to multigrid is given in [37]. Model codes and a detailed description of application to (elasto-)hydrodynamic and non-adhesive contact problems in the classic boundary integral formulation are given in [19]. A wide range of problems is treated in [38]. For completeness, in Appendix C a basic description of geometric multigrid is given with implementation details for the Navier–Cauchy equations with the contact problem boundary condition as used in this research.

4. Verification and model results Two model problems are considered: firstly, the (classic) academic case of contact between a flat homogeneous material and a spherical (paraboloid) indenter, and secondly, the contact problem between a spherical indenter and a homogeneous material with soft and hard coating respectively.

4.1. Homogeneous material The analytic solution for a semi-infinite homogeneous bulk in contact with a paraboloid indenter of radius of curvature R is the Hertzian solution [12] of a semi-ellipsoid contact pressure with maximum: ph ¼

3F ; 2π a2h

ð10Þ

Fig. 1. Flow diagram of a Full MultiGrid (FMG) algorithm starting at a coarse grid passing through finer grids (levels) to a target grid at each level using the coarser grids for convergence acceleration and eventually arriving at a solution at the target grid. Shown is the case with 4 gridlevels, and a single coarse grid correction cycle type V per level. The circles indicate single grid iterations to smooth the error. The arrows up refinement/correction from coarse to fine mesh. The arrows down coarsening of a smooth error. The double circles are the solutions achieved at different grids, see Appendix C.

236

and a contact radius
1=3 3FRð1  ν2s Þ ; ah ¼ 4Es

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

ð11Þ

where νs and Es are Poisson's ratio and Young's modulus of the homogeneous bulk respectively. F is the external load. The equations were non-dimensionalized using the Hertzian solution, i.e. the x; y; and z coordinates scaled on ah, and the pressure and stresses on the maximum Hertzian contact pressure. The non-dimensional problem has been solved using different sizes of the 3D domain: ½6:0; 6:0; 3:0nah , ½6:0; 6:0; 6:0nah and ½8:0; 8:0; 8:0nah . Young's modulus is chosen as Es ¼ π =2nð1  ν2 Þ. For each case the problem is solved on grids with decreasing mesh size ranging from 1/2 to 1/128 on the finest grid. The coarsest grid in the multigrid algorithm in each case had a mesh size of unity, and the number of gridlevels ranged from 2 to 8. In Fig. 2 the relative difference Θ between the maximum pressure P0 at the surface from the numerical solution of the 3D contact problem on each of the domains, and the maximum Hertzian value, is shown as a function of the number of nodes per unit length. The result shows that with decreasing mesh size the difference reduces to a fixed value which indicates the (small) difference between the solution for the 3D finite domain considered here, and the Hertzian solution which holds for an infinite domain. For the smallest domain the difference is about one percent. As can also be seen, increasing the domain decreases the difference between the converged solution and the semi-infinite body solution. The results for each domain show the asymptotics of grid convergence with a slope of almost 2 (comparable to the curve of equation y ¼ x  2 ) which is consistent with the second order discretization used. In fact, owing to the effect of the singularity in Δp in the Hertzian solution the slope should be between 1.5 and 2. To conclude, this study shows that the finite domain solution converges to the infinite halfspace solution with increasing domain size as should be expected.

4.2. Coated materials For a detailed discussion on the effect of coatings the reader is referred to [39]. The contact problem considered is depicted in Fig. 3: a spherical indenter of radius R in contact with a solid which consists of a coating of thickness t on top of a material with homogeneous properties: elastic modulus Es and Poisson ratio νs . The elastic modulus of the coating is Ec and its Poisson ratio νc . Stevanovic et al. [40] presented a (curve-fit) solution for the

Fig. 3. Elastic contact of a sphere with a layered bulk.

contact radius, based on the work of Chen and Engel [41]:
pffiffiffiffiπ =4 !! 1=4 t α a ¼ aS þ ðaC  aS Þ 1  exp  π ; a0 with a0 ¼ aS þ ðaC  aS Þ 1  exp  π

pffiffiffiffi π =4 !! 2t α : aS þ aC

ð12Þ


1=4

ð13Þ

where aC and aS are the contact radius for the case of an infinitely thick ðt-1Þ and zero (t ¼0) coating thickness respectively. The parameter α is defined as the ratio

α

aC aS

ð14Þ

Eq. (12) is supposed to be accurate for α o 5. The cases considered in this section all are in this regime. The domain size used for the computation is ½8:0; 8:0; 4:0nah or ½8:0; 8:0; 6:0nah (thick coatings), so that the boundary conditions do not affect the solution even for the case of the entire volume consisting of the softest material used for the coating. The ratios of Young's modulus considered are Ec =Es ¼ 0:05; 0:1; 0:25; 0:5; 1:0; 2:0 and 3.0 for the case of (dimensionless) coating thickness values of t=ah ¼ 0:25; 0:5; 1:0 and 2.0. The value of Poisson's ratio is ν ¼ νc ¼ 0:3 for all materials. The externally applied dimensionless load F ¼ 2nπ =3, so that for Ec =Es ¼ 1 the problem is the equivalent of the Hertzian problem. Fig. 4 shows the contact pressure profiles solved for different cases. Comparing the values of the contact radius from the computed solutions with Eq. (12) from [40], the maximum difference is 1.3% for Ec =Es ¼ 0:05, which shows the curve fit formula to be quite accurate. Boffy et al. [27] used a multigrid method to solve the elastic problem combined with a conjugate gradient algorithm to determine the contact pressure for the same problem. The present results are also in very good agreement with those findings, as should obviously be the case. However, in the present approach the solution is obtained more efficiently as the entire problem is solved in an amount of work that is the equivalent of only one or a few times solving the elastic problem itself.

5. Heterogeneous material applications In this section results are presented for more demanding cases representative for the applications for which the algorithm has been developed. First the case of a homogeneous material with a cluster of inclusions is considered. 5.1. Cluster of inclusions Fig. 2. Relative error Θ between the numerical solution P0 and the analytical one Ph. The number of level used for the simulations varies between 2 and 8.

Consider a homogeneous material with a local cluster of inclusions in a rolling contact configuration. When located sufficiently deep

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

237

compared to the case of a single inclusion having the same mechanical properties (Ei, νi ). Fig. 5 shows the topology cross section in the central plane y=ah ¼ 0. In 3D a cluster of 20 spherical inclusions is arranged on a sphere of radius R=ah ¼ 0:125. The radius of the inclusions composing the cluster is r i =ah ¼ 0:04. Results for the cluster identified with C will be compared with results for a single large inclusion (identified with SLI) with its center located at the same location as the center of the cluster. The location of the center of the cluster will be varied in depth with d ¼ z=ah . Finally a single small inclusion is considered. This case, identified as SSI, represents the case where top inclusion of the cluster is considered, which is easily generated from the cluster by setting Young's modulus of the other 19 small inclusions at Ei =E ¼ 1. The domain size used for the computation is ½6:0; 6:0; 3:0nah . The calculation was done with a mesh size h ¼ 1=128 which gives about 225  106 grid points. Each inclusion is meshed with more than 500 grid points to ensure the accuracy of the solution.

below the surface this cluster will not significantly affect the fatigue life of the contact, as it does not alter the pressure distribution. It may still affect the maximum von Mises stress in a local region but the overall maximum may not increase above the standard maximum for a Hertzian pressure. However, when located closer to the surface the pressure and the stress distribution may significantly be altered. To develop advanced failure criteria, allowable size and distribution of pores in ceramics, understanding of the effect of inclusions and clusters of inclusions in the actual contact configuration is important. The effect on the stresses, i.e. the fatigue life, will depend on the properties of the inclusions in the cluster relative to the bulk, i.e. being harder or softer than the bulk material, on geometry, and on the location. For engineering practice one may wonder if certain configurations of clusters may be viewed as a single “equivalent” inclusion. In this section some results are presented of the contact pressure and the sub-surface stress field for cases of a contact problem with a model cluster of inclusions in the subsurface material. The results are

2

2 E /E=0.05

E /E=0.05

E /E=0.1

E /E=0.1

E /E=0.25

E /E=0.25

c

c

c

1.5

c

1.5

c

E /E=0.5

c

E /E=0.5 c

P/P0

P/P0

c

Ec/E=1.0 E /E=2.0 c

1

E /E=3.0 c

0.5

0 −3

E /E=1.0 c

Ec/E=2.0

1

E /E=3.0 c

0.5

−2

−1

0

1

2

3

x/a

0 −3

−2

−1

0

1

2

3

x/a

0

0

2 E /E=0.05 c

E /E=0.1 c

E /E=0.25

1.5

c

E /E=0.5

P/P0

c

Ec/E=1.0 E /E=2.0 c

1

E /E=3.0 c

0.5

0

−3

−2

−1

0

1

2

3

x/a

0

Fig. 4. Contact pressure in the presence of a coating having different thicknesses t. From left to right, top to bottom t=ah ¼ 0:25; 0:5; 1:0; 2:0.

Fig. 5. Elastic contact of a sphere with a homogeneous bulk in the presence of single large inclusion (SLI), a cluster of inclusions (C) or a single small inclusion (SSI).

238

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

cluster is almost the same for z=ah small. However, the total width of the region in which the contact pressure is affected for the cluster is comparable to the size of the disturbance for a large single inclusion. With increasing depth, the effect of the cluster reduces faster than the effect of the single large inclusion. This is explained by the smaller ratio of inclusion radius to depth. For example at a depth z=ah ¼ 0:4 the pressure reduction induced by a single large inclusion is still about 10%, whereas that of the cluster only about 2%. For z=ah 40:5 the effect of the cluster has almost completely disappeared. In Table 1 also the values of the maximum von Mises stress below the surface are shown for the different cases. The von Mises stress field in the central plane for the case of a cluster, a single small inclusion, and a single large inclusion, at a depth z=ah ¼ 0:25, are shown in Fig. 7. The graphs illustrate where, for these cases, the maximum von Mises stress occurs, i.e. on each side of the (single) small and large inclusion, where the inclusions themselves are low stress regions, and, for this depth, at the location of the inclusions below the top for the cluster. In the case of a cluster, the density of inclusions maintains the maximum level of stress almost constant in the subsurface independent of the depth. This maximum is indeed carried successively by the different inclusions composing the defect. Comparison with the single top inclusion of this cluster shows Table 1 Contact pressure Pð0Þ=P h and maximum von Mises stress σ VM =P h as a function of the cluster (inclusion) depth for Ei =E ¼ 0:4. The center of the contact is at ðx=ah ¼ 0:0; y=ah ¼ 0; z=ah ¼ 0Þ. Depth z=ah

Pð0Þ=P h C–SSI–SLI

Max ðσ VM =P h Þ C–SSI–SLI

0.2 0.25 0.4 0.55 1

0.761–0.818–0.570 0.910–0.950–0.652 0.984–0.996–0.885 1.000–1.000–0.956 1.000

0.803–0.662–0.810 0.823–0.659–0.860 0.825–0.757–0.930 0.835–0.797–0.911 0.621

Fig. 6. Solution of the contact problem in the presence of a soft C (top), SSI (middle), SLI (bottom). z=ah represents the depth of the center of both the cluster and the SLI.

The contact problem is solved for F ¼ 2nπ =3 which for Ei =E ¼ 1 is identical to the homogeneous Hertzian problem. The center of the indenter is located at ðx=ah ¼ 0; y=ah ¼ 0; z=ah ¼ 0Þ. Both pressure and stress profiles are investigated for soft ðEi =E ¼ 0:4Þ and hard ðEi =E ¼ 2:5Þ clusters of inclusions (C) and inclusions (SSI and SLI). The value of Poisson's ratio is ν ¼ νi ¼ 0:3 for all materials.

5.1.1. Soft inclusions Fig. 6 shows centerline pressure profiles along x=ah for the cases when the center of the cluster/single large inclusion is located at depths z=ah ¼ 0:2; 0:25; 0:4; 0:55, and 1, which is the homogeneous case. The calculated pressure on the center of the contact and the maximum von Mises stress in the bulk for these configurations are given in Table 1. These results show that the different configurations of soft inclusions all have a reduced maximum contact pressure, as expected. Comparing the results for a single large and single small inclusion with the result for the cluster shows that the effect on the contact pressure is dominated by the top inclusion in the cluster. The pressure drop for a single small inclusion and for a

Fig. 7. Von Mises stress in the presence of a soft C (top), SSI (middle), SLI (bottom) located at z=ah ¼ 0:25.

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

239

that this inclusion is carrying the maximum when the cluster is located deep in the subsurface ðz=ah ¼ 0:55Þ. In the case of a single large inclusion it can be observed that the maximum value of the von Mises stress is reached for a depth close to z=ah ¼ 0:4. The stress intensity is higher compared to the one obtained with the cluster (8%) which means that this configuration is more critical in terms of fatigue life. The figure shows that the topology of the von Mises stress field differs for the different cases. The cluster represents a region of locally relatively high stresses, in particular at the inclusions just next to the centerline for the cluster on the surface and on the deep side. The top inclusion itself is a region of small stress. For the small inclusion the inclusion itself is a low stress region, but at its side two (local) maxima occur. The maximum value here is still to the sides of the single small inclusion, and to the sides of the single large inclusion.

5.1.2. Hard inclusions Next results are presented for the case the inclusions are harder than the bulk material. Table 2 shows the pressure at the center of the contact and von Mises stress for the different cases. The results show that when hard defects are deep below to the surface, the cluster behaves almost exactly like the single top inclusion located on its top. When the cluster is shallow, the intensity of the von Mises stress is close to the one of a single large inclusion. However, the maximum pressures in these cases are not comparable as far as the overstress created by a shallow large inclusion would lead to a premature failure in the material. As a conclusion to this study, it has been observed that a cluster of soft inclusion could be modeled as a single soft large inclusion when it is located shallow to the surface and that a cluster of hard inclusion could be modeled as a single hard small inclusion when it is located deeper. This result may help in the design of new materials.

Fig. 8. Granular matter with an ellipsoidal inclusion and fractured surface of silicon nitride Si3N4.

5.2. Polycrystalline material Finally, an application representative for ceramic materials is presented. The subsurface topology consists of homogeneous grains with an ellipsoidal inclusion as an approximation of a needle. The grains are separated by a layer. The results are presented to show the ability of the method to deal simultaneously with both irregular geometries like grains and property variations. Fig. 8 shows the modeling of the problem and an illustration of a fractured surface of ceramic (silicon nitride Si3N4). 1960 grains are designed in a ½6:0; 6:0; 3:0nah bulk using a 3D Voronoi tessellation. The inclusion (needle) center is at ðx=ah ¼ 0; y=ah ¼ 0; z=ah ¼ 0:5Þ. The domain is covered with a uniform grid of some 225 millions points, allowing each grain to contain in excess of than 105 points. Between the grains a layer of glue is modeled which in total represents 10% of the domain volume. Table 2 Evolution of the contact pressure Pð0Þ=P h and maximum von Mises stress σ VM =P h as a function of the defect depth for hard defects. The center of the indenter is located at ðx=ah ¼ 0; y=ah ¼ 0; z=ah ¼ 0Þ. Defect depth z=ah

Pð0Þ=P h C–SSI–SLI

Max ðσ VM =P h Þ C–SSI–SLI

0.2 0.25 0.4 0.55 1

1.353–1.239–1.917 1.123–1.070–1.456 1.041–1.021–1.287 1.024–1.014–1.194 1.000

1.336–0.937–1.358 1.345–1.142–1.314 1.359–1.327–1.304 1.353–1.355–1.274 0.621

Fig. 9. From top to bottom: solution of the contact problem for a polycrystalline material with an ellipsoidal inclusion – case 1, central indenter, Eg ¼ 2=5E – case 2, non-central indentor ðx=ah ¼ 0; y=ah ¼ 0:5Þ, Eg ¼ 2=5E – case 3, non-central indenter, softer glue Eg ¼ E=5 and top glue layer.

240

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

Fig. 10. From top to bottom: von Mises stress in the central planes (x¼ 0 and y¼ 0) for the contact problem of a polycrystalline material with ellipsoidal inclusion – case 1, central indenter, Eg ¼ 2=5E – case 2, non-central indenter ðx=ah ¼ 0; y=ah ¼ 0:5Þ, Eg ¼ 2=5E – case 3, non-central indenter, softer glue Eg ¼ E=5 and top glue layer.

5.2.1. Case 1 In this case we consider a volume which is entirely composed of grains and glue. Young's moduli of the glue and the inclusion are chosen 2/5 times and 5 times higher than the one of the grains respectively. Poisson's ratio is taken the same for the different phases, i.e. ν ¼ 0:3. Fig. 9 (top) shows the pressure profile at the top surface. As expected the grain and glue phases create discontinuities whose shape and intensity depend on the microstructure. Fig. 10 (top) shows the von Mises stress in both x=ah ¼ 0 and y=ah ¼ 0 planes. The maximum stresses occur at the rim of the inclusion, where the different phases are located. The position of the grains determines the stress path along the inclusion. The interstitial glue forms low stress (relaxed) regions and consequently the glue effectively forms stress paths barriers between the grains. 5.2.2. Case 2 This case is the same as the previous one excepted that the indenter is located at ðx=ah ¼ 0; y=ah ¼ 0:5; z=ah ¼ 0Þ. Fig. 9 (middle) shows the pressure profile solved in this configuration. Fig. 10 (middle) shows the von Mises stress in both x=ah ¼ 0 and y=ah ¼ 0 planes. The surface pressure clearly reflects the results of the shift of the indenter center. The maximum stresses occur at the same location as for case 1, but are at a lower level. The overall view of the stress field is quite similar. The differences are in the details of the grains and inclusion. 5.2.3. Case 3 This case considers the same volume, and location of the indenter as in case 2, but now a layer of glue is modelled between

the grains and the top surface. This layer will play the role of a soft coating between the surface and the microstructure. Young's modulus of the glue is now chosen lower, i.e. 1/5 times the value of the grain material. The indenter is located at ðx=ah ¼ 0; y=ah ¼ 0:5; z=ah ¼ 0Þ as for case 2. Fig. 9 (bottom) shows the surface pressure profile solved for this case. Fig. 10 (bottom) shows the von Mises stress in both x=ah ¼ 0 and y=ah ¼ 0 planes. The surface pressure distribution shows the surface glue layer mollifies the variations and “absorbs” high frequency variations resulting from the deeper granular structural variation. The results show that the maximum stresses occur at the interface between the inclusion, grains and glue and reach a higher intensity. This difference is due to Young's modulus of the glue which has been chosen lower in this configuration. In the layer of glue the level of stress reminds low. However, as for a single inclusion, the most critical load position is not often located straightly upper this one but at a certain distance that depends on the defect properties. This kind of studies, combined with fatigue criterion, may help in the understanding of the failure in the case of advanced complex materials.

6. Conclusion An efficient geometric multigrid algorithm has been developed to solve 3D elastic contact problems for heterogeneous materials with realistic subsurface topology. The examples presented in this paper illustrate the potential of the method to handle (local) property variations in value as well as in an (irregular) microstructure. The flexibility to use dense grids used allows to precisely

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

describe the pressure profiles and give an accurate solution for the stress fields. The algorithm is well suited to be used for contact problem analysis in for a wide variety of applications ranging from classic rolling contacts to tire road contacts and biomedical applications. It can greatly contribute to the development of new (rolling contact) fatigue models and aid the development of design specifications composite and polycrystalline materials in contact applications. In this paper a rigid indenter body was assumed. However, this is not a limitation of the approach. One can extend the presented algorithm to the case of two heterogeneous materials with different topologies in contact. Finally, the method can be generalized and extended to elasto-hydrodynamic lubrication contacts as well by introducing the Reynolds equation as the boundary conditions to determine the pressure and film thickness profiles.

where uh, vh and wh indicate the exact solution of the approximate system of equations on the grid with meshsize h. In this study the 2 following Oðh Þ approximations were used for each of the terms at the gridpoint with index ði; j; lÞ: Lhx;1 ¼

Q i þ 1=2;j;l uhiþ 1;j;l  ðQ i þ 1=2;j;l þ Q i  1=2;j;l Þuhi;j;l þQ i  1=2;j;l uhi 1;j;l 2

Lhx;2 ¼ Lhx;3 ¼ Lhx;4 ¼ Lhx;5 ¼

Acknowledgments

Lhx;6 ¼

The authors gratefully acknowledge the support of SKF ERC, Nieuwegein, The Netherlands. The project is also partly carried out in the framework of the innovation program GO Gebundelde Innovatiekracht, and funded by the European Regional Development Fund, Regio Twente and Provincie Overijssel. The project partners Apollo Tyres Global R&D, University of Twente (Tire-Road Consortium), Reef Infra, Stemmer Imaging and the Provincie Gelderland are gratefully acknowledged.

241

Lhx;7 ¼

h

λ

λ

h h h h i þ 1;j;l ðvi þ 1;j þ 1;l  vi þ 1;j  1;l Þ  i  1;j;l ðvi  1;j þ 1;l  vi  1;j  1;l Þ 2

4h λi þ 1;j;l ðwhiþ 1;j;l þ 1 whiþ 1;j;l  1 Þ  λi  1;j;l ðwhi 1;j;l þ 1  whi 1;j;l  1 Þ 2

4h

μi;j þ 1=2;l uhi;j þ 1;l  ðμi;j þ 1=2;l þ μi;j  1=2;l Þuhi;j;l þ μi;j  1=2;l uhi;j  1;l 2

h

μi;j;l þ 1=2 uhi;j;l þ 1  ðμi;j;l þ 1=2 þ μi;j;l  1=2 Þuhi;j;l þ μi;j;l  1=2 uhi;j;l  1 2

h

μi;j þ 1;l ðvhiþ 1;j þ 1;l  vhi 1;j þ 1;l Þ  μi;j  1;l ðvhiþ 1;j  1;l  vhi 1;j  1;l Þ 2

4h

μ

μ

h h h h i;j;l þ 1 ðwi þ 1;j;l þ 1  wi  1;j;l þ 1 Þ  i;j;l  1 ðwi þ 1;j;l  1  wi  1;j;l  1 Þ 2

4h

where the “half” indexes indicate that values are taken at locations in the center between two gridpoints. The discrete equations for the y and z directions in the interior points can be derived in exactly the same way.

Appendix A. Discrete equations

A.2. Boundary conditions

All equations were discretized with second order accuracy on a uniform grid using a collocated approach, i.e. the equations and unknowns are defined at each point of the grid, see [27,35]. Higher order discretizations can be developed using, e.g., the so-called compact schemes, see [42]. When the material is near incompressible ðν-0:5Þ a different approach is required. In that case the discretization needs to be done using a staggered grid to avoid spurious oscillations as for the Stokes equations describing incompressible slow viscous flow, e.g., see [43].

The boundary conditions for the plane where the contact occurs ðz ¼ 0Þ are no stress in the x and y directions, and when there is contact ðf ðx; yÞ ¼ wðx; y; 0ÞÞ the normal stress σ zz balances the contact pressure pðx; yÞ. This leads to the following equations:

A.1. Interior equations The Navier–Cauchy equations have been described in Section 2. Defining Q ¼ λ þ2μ the equation for the x direction is



 ∂ ∂u ∂ ∂v ∂ ∂w ∂ ∂u Q þ λ þ λ þ μ ∂x ∂x ∂x ∂y ∂x ∂z ∂y ∂y |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} 1 2 3


 4 ∂ ∂u ∂ ∂v ∂ ∂w þ μ þ μ þ μ ¼0 ð15Þ ∂z ∂z ∂y ∂x ∂z ∂x |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} 5

6

7

which can be written as Lx ðu; v; wÞ ¼

7 X

Lx;s ðu; v; wÞ ¼ f x

ð16Þ

s¼1

with Lx;s ðu; v; wÞ being the operators indicated in (15) and fx the right-hand side of the x equation. Although this right-hand side is zero here, the description is given for a general right-hand side, which is relevant in view of the application in a multigrid algorithm. On a uniform grid with meshsize h in each of the dimensions, a discrete approximation to Eq. (16) for each interior gridpoint with index ði; j; lÞ can be written as Lhx ðuh ; vh ; wh Þijl

¼

7 X s¼1

Lhx;s ðuh ; vh ; wh Þijl

¼

h f x;ði;j;lÞ

ð17Þ

∂u ∂w þ ¼ 0; ðaÞ ∂z ∂x ∂v ∂w þ ¼ 0; ðbÞ ∂z ∂y f ðx; yÞ  w ¼ 0; ðcÞ
 ∂w ∂u ∂v þ þ p ¼ 0; ðλ þ 2μÞ þ λ ∂z ∂x ∂y

ðdÞ

ð18Þ

where the third equation cancels when the gap is open. In that case the third and fourth equations are replaced by the fourth equation with p ¼ 0, hence, the problem is a “complementarity” problem. The appropriate formulation of the (normal) contact condition is p Z 0, f ðx; yÞ  w Z0, and pðf ðx; yÞ  wÞ ¼ 0 (Hertz– Signorini–Moreau). Although in this paper we only consider normal contact pressure, including a friction boundary condition in these equations can be done easily. Let the plane z¼0 have grid index 0 in the z direction, and the gridplanes below it indexes 1 and 2. A second order approximation to the equations in the gridpoint with indexes ði; j; 0Þ of the top surface is 1 1 ð1:5uhi;j;0  2uhi;j;1 þ 0:5uhi;j;2 Þ þ ðwhiþ 1;j;0  whi 1;j;0 Þ ¼ 0 ðaÞ h 2h 1 1 ð1:5vhi;j;0  2vhi;j;1 þ 0:5vi;j;2 Þ þ ðwhi;j þ 1;0 whi;j  1;0 Þ ¼ 0 ðbÞ h 2h f i;j;0  whi;j;0 ¼ 0 ðcÞ 1 ðλi;j;0 þ2μi;j;0 Þ ð1:5whi;j;0  2whi;j;1 þ 0:5whi;j;2 ÞÞ h þ

λi;j;0 2h

ðuhiþ 1;j;0  uhi 1;j;0 Þ þ

λi;j;0 2h

ðvhi;j þ 1;0  vhi;j  1;0 Þ þ phi;j;0 ¼ 0

ðdÞ ð19Þ

This ð4  4Þ system holds only for otherwise the gap is open, the third equation cancels, leaving a ð3  3Þ system and phi;j;0 ¼ 0. This phi;j;0 4 0,

242

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

is to obtain a solution that satisfies the Hertz–Signorini–Moreau condition. Regarding the other boundary planes of the domain, the boundary conditions for the bottom plane u ¼ v ¼ w ¼ 0 are trivial. The boundary conditions for the other faces of the domain can be derived in the same way as above without the equations determining the contact. This results in a ð3  3Þ system of equations. The details are left to the reader. The last equation to be considered is the force balance equation. Assuming a piecewise constant pressure for each point ði; j; 0Þ in the contact boundary plane yields a second order approximation to Eq. (9): XX 2 h phi;j;0 ¼ f i

ð20Þ

j

Appendix B. Relaxation In this appendix a Gauss–Seidel relaxation is described for the ~ h denote the current system of equations. Let u~ h , v~ h and w approximation to uh, vh, and wh. The specific value at a gridpoint is denoted by the symbols with a subscript ði; j; lÞ.

B.1. Interior equations A system of equations is characterized by its determinant scalar operator. The determinant of the Navier–Cauchy equations is the biharmonic and therefore the system is elliptic. When designing a relaxation for a system of equations one needs to take into account the coupling between the equations. When λ and μ are such that the problem is sufficiently far from the incompressible limit ðν ¼ 0:5Þ the equations for the different directions are weakly coupled, meaning that the terms related to v and w in the x equation are of secondary importance relative to the terms related to u. This implies that the x equation can be considered as an equation for u, the y equation as an equation for v, and the z equation as an equation for w, and one can relax the equations consecutively, i.e. one after the other. Below the relaxation for the x equation is described. A relaxation sweep for this equation consists of scanning the grid in a prescribed (lexicographic) order, and at each point ði; j; lÞ the current approximation to the unknown u~ hi;j;l is updated: u i;j;l ¼ u~ i;j;l þ r hx

∂Lhx ∂uhi;j;l

!1 ð21Þ

rhx

where is the dynamic residual (in the point ði; j; lÞ) before making the change, defined as h h ~ hÞ r hx ¼ f x  Lhx ðu^ ; v~ h ; w

ð22Þ h

The residual is computed using u^ which is equal to u~ for all points which have not yet been relaxed, and to u h for all points which have already been changed. So, as is characteristic for Gauss–Seidel relaxation, changes already made are taken into account. For a linear problem Eq. (21) simply states that a change is made to the current approximation to uh in the gridpoint such that the discrete equation at that point, fixing the values in the surrounding points, is exactly satisfied after the change is made. For a non-linear equation it is a single Newton iteration step to solve the equation at that gridpoint by changing the value of the unknown at that point only. From Eq. (17) and the definitions of h

the operators Lhx;s it can be seen that ∂Lhx ∂uhi;j;l

!

¼

1 2

h

 ðQ i þ 1=2;j;l þ Q i  1=2;j;l þ μi;j þ 1=2;l þ μi;j  1=2;l þ μi;j;l þ 1=2 þ μi;j;l  1=2 Þ



ð23Þ B.2. Boundary equations For the sides of the domain with stress free boundary condition the relaxation in the boundary plane can be carried out in the same way as for the interior, the boundary plane is scanned point by point and at each point the boundary equation associated with u is solved for u, the one for v is solved for v, etc. The detailed equations can be derived in the same way as was done for the interior above, except that now the discretized boundary operators should be taken. The equations can be relaxed consecutively because the coupling between the equations is weak, i.e. in the equation for u the terms associated with v and w are not dominant, so that when making changes such that one equation is solved the other equations at the same point and its neighbors are not disturbed too much. A formal analysis can be carried out using the so-called local mode analysis, see [19]. This is not the case for the boundary plane in which contact occurs. Here the complementarity and associate change of the number of boundary equations creates a strong coupling between the equations at the same point. In this case a simple consecutive relaxation will diverge. The solution is to simultaneously solve for all equations at the boundary point. The resulting relaxation is referred to as coupled or collective relaxation. This implies that a small 4  4 system of equations needs to be solved at each gridpoint in the boundary plane. For the point ði; j; 0Þ the system of equations expressed in terms of the corrections to be applied can be written as 0 1 1 0 1 0 r h1 a1;1 a1;2 a1;3 a1;4 δu B hC C B C B Br C B a2;1 a2;2 a2;3 a2;4 C B δv C B 2C C B C B ð24Þ C B a3;1 a3;2 a3;3 a3;4 C B δw C ¼ B B r h3 C A @ A @ @ A a4;1 a4;2 a4;3 a4;4 δp i;j r h4 i;j i;j

and r the residual of equation l ðl ¼ 1; 2; 3; 4Þ of where δ the system (24) corresponding to the node ði; jÞ in the boundary plane. The coefficients am;n represent the weights of unknown n in the m equation. The system can be solved directly in a small number of operations. Indeed, this is a computationally more expensive approach than the consecutive relaxation. However, as this is done only for the boundary plane the effect on the total operation count of the relaxation is small. Subsequently, the computed changes are applied: u ¼ u hi;j  u~ hi;j

l

u hi;j ¼ u~ hi;j þ ωδu i;j where ω is an underrelaxation parameter, and in the same way new approximations v h , w hi;j , and p hi;j are computed. After the changes are applied the sign of the pressure is verified. When the new pressure value is negative, indicating the surfaces at this point are not in contact, the pressure is set to zero, the third equation is removed and changes are applied according to the solution of the resulting 3  3 system for a stress free surface. From the discrete equation the system of equations can easily be derived. For example, the full 4  4 system is given by 0

3  2h

B B 0 B B 0 @ 0

0

0

3  2h 0

0 1

0

3  2h

1 0

0

1m

1 δr 1i;j δui;j m þ 1 B 2 C C B δv C B δr i;j C 0 C B i;j C B C CnB C ¼B 3 C B δwi;j C B δr i;j C 0C A @ A @ A δpi;j 1 δr 4

0

i;j

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

B.3. Relaxation: force balance equation The force balance condition, Eq. (9), determines the value of the separation constant δ in (24). This equation is of the type global constraint because it does not just determine the solution in a single or few points but has an effect on the entire solution. A complication is that the value of δ does not appear in the equation itself. This can be resolved by an iteration as follows, see also [19]: 0 1 X 2 h ð25Þ δnew ¼ δold þ ω@f  h pi;j A i;j

where ω is a relaxation parameter. When using a single grid this change to δ can be made every so-many relaxations, where for stability it is needed to keep the changes limited. In a multigrid algorithm a global constraint can be dealt with in two ways. This is explained in Appendix C.

order of a few mesh sizes, i.e. high frequency error components. As a result, after a few iterations the error v~ h is smooth relative to the mesh size. Rather than continuing the iterative process on grid h 2 with a slow convergence rate 1  Oðh Þ in a multigrid algorithm the smoothness of the error is exploited. Let rh denote the residual defined as h

r h ¼ f  Lh u~ h :

~ in all points of the domain, and p~ in the Let u~ , v~ , and w contact plane, be the current approximation to the solution: a Gauss–Seidel step now consists of h

h

 Collective relaxation of all points in the boundary plane z ¼0 to obtain u h , v h , w h , and p h in these points.

 Sequential relaxation of the boundary planes with zero stress 

h h

By definition f ¼ L u so this equation can be written as r h ¼ Lh uh  Lh u~ h ;

ð29Þ

which for the case of a linear operator can be written as Lh vh ¼ r h

ð30Þ

As the error is smooth this equation can, with little loss of information, be represented on a coarser grid, say with meshsize H ¼ 2h:

condition on the sides, to obtain u h , v h and w h in the points of these planes. Sequential relaxation of the interior in lexicographic order: ○ Compute u h in all interior points ði; j; kÞ from the equation for the x direction, ○ Compute v h in all interior points ði; j; kÞ from the equation for the y direction, ○ Compute w h in all interior points ði; j; kÞ from the equation for the z direction.

ð31Þ

where L is a suitable approximation to L defined on grid H and IH h is a restriction operator also to be defined. Often LH is simply the discretization of the continuous operator L on the grid H, and IH h the transpose of a linear interpolation in each dimension, which gives a weighted averaging. Next, Eq. (31) is solved for vH on the coarse grid. As the equation is essentially the same as the fine grid equation, except for the fact that it is defined on a coarser grid, this can be done using the same iterative scheme. The advantage is that it will converge faster, at a rate 1  OðH 2 Þ and the number of nodes is also smaller by a factor ðH=hÞd when d is the dimension of the problem. Hence, solving the equation can be done much more efficiently. Once vH (or a sufficiently accurate approximation to it) has been solved it can be used to correct the fine grid approximation according to H

h

ð28Þ h

h LH vH ¼ I H hr

B.4. Summary of relaxation h

243

h

u h ¼ u~ h þ I hH vH IhH

Carrying out these steps once will be referred to as a single grid iteration or relaxation sweep. This process, alternated once every so many steps with a relaxation of the force balance equation according to (9) by itself, completes a numerical solution method. It is stable but it is has a very slow convergence speed owing to the fact that it is only efficient in reducing errors with a wavelength comparable to the 2 mesh size. Smooth components are only reduced at a factor 1  Oðh Þ per sweep leading to very slow asymptotic convergence. This slowness is overcome by using multigrid techniques.

Appendix C. Multigrid: formulation and implementation In this appendix a generic description of a multigrid algorithm in symbols and equations is given, see also [44,19]. The problem to be solved is written as L h uh ¼ f

h

ð26Þ

where uh stands for the solution to be determined on a grid with meshsize h. Lh stands for the discrete operator, which may also be non-linear, i.e. depend on uh. So, the notation Lh uh should not be read as a matrix–vector multiplication but rather as Lh applied to uh. Let u~ h stand for the approximation to uh obtained as a result of a number of iterations using a process such as Gauss–Seidel relaxation. The error in the approximation is given by v~ h ¼ uh  u~ h

where stands for an interpolation from the coarse grid to the fine grid. As the correction will be smooth generally a linear interpolation in each dimension is sufficiently accurate. When the coarse grid is still too fine for the relaxation to be an efficient solver the idea can be applied recursively, i.e after a few relaxations on the coarse grid, one transfers the problem to an even coarser grid for the correction of the solution of the remaining smooth error. This process is repeated until a grid is reached at which the problem can be solved using only a few iterations. Subsequently the coarse grid solution of the error is used to correct the next finer grid, and so on until the target grid is reached again. This process is called a coarse grid correction cycle. For the case of 4 grids (levels) the flow diagram of a cycle is shown in Fig. 11. Owing to the shape of the cycle it is referred to as a Vðν1 ; ν2 Þ cycle, where ν1 and ν2 are the number of relaxations carried out on each grid before coarsening and after returning to the grid respectively. The potential error reduction of a coarse grid correction cycle is μ ν1 þ ν2 where μ is the asymptotic smoothing factor defined as the largest amplification factor of components which can only be seen on the fine grid. For elliptic problems and properly designed

ð27Þ

Assume that, as is often the case, the relaxation processes are efficient only in reducing components with a wavelength of the

ð32Þ

Fig. 11. Flow diagram of a V ðν1 ; ν2 Þ cycle for the case of 4 grids (levels).

244

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

relaxation typically μ ¼ 0:5  0:8. As a result a cycle with only a few relaxations per grid has a grid independent convergence factor of an order of magnitude, and only a few cycles are needed to obtain a converged solution. Neglecting the work done in the intergrid transfers the computational work of a cycle is approximately W

ðν1 þ ν2 Þ 1 ðH=hÞ  d

ð33Þ

Wu

where Wu is the amount of work of a single relaxation on the finest grid. For a 3D problem with H=h ¼ 2 this is only 8/7 times the work of the ν1 þ ν2 relaxations actually carried out on the finest grid. Hence, the cycle represents a huge increase in efficiency compared to single grid relaxations and when combined with using the coarser grids for generating an accurate first approximation, see Fig. 11. It is beyond the scope of this section to discuss all aspects of efficiency and algorithm design. Further details can be found in the vast multigrid literature. When the problem is non-linear the formulation needs to be changed. In that case (30) does not hold. However, using (27) and (29) it follows that Lh ðu~ h þ vh Þ ¼ r h þ Lh u~ h :

ð34Þ

This equation is, by definition, exact for any problem and reduces to Eq. (30) when the operator is linear. The equation contains as unknown now sum of the error and the current approximation. A coarse grid representation of (34) is h ~h H h LH u^ ¼ LH I H h u þ Ih r ;

ð35Þ

H where the coarse grid variable u^ is defined as H u^ ¼ I hH u~ h þ vH :

ð36Þ

As Eq. (35) is the same as the fine problem, so it can be solved using the same iterative procedure. Only the right-hand side is different, the meaning of the variable solved is different, and it H contains less unknowns. Once u^ has been solved on the fine grid approximation is corrected according to H

~h u h ¼ u~ h þ I hH vH ¼ u~ h þ I hH ðu^  I H h u Þ:

ð37Þ

This scheme for the non-linear problem is referred to as the Full Approximation Scheme (FAS), see [44]. The formulation for the linear problem is called the Correction Scheme. For more details, including the implementation in model codes the reader is referred to, e.g., [19,38,45] and many other books and papers in the multigrid literature. C.1. Implementation for Navier–Cauchy equations To conclude this section some specific details regarding the implementation of multigrid techniques for the Navier–Cauchy equations are given. C.1.1. Relaxation As explained in Appendix B, a system of equations is characterized by its determinant scalar operator. The determinant of the Navier–Cauchy equations is the biharmonic and therefore the system is elliptic. When designing a relaxation for a system of equations one needs to take into account the coupling between the equations, see [45,38]. For the Navier–Cauchy equations when the values of the Lamé parameters are sufficiently far from the incompressible limit simple one point Gauss–Seidel relaxation applied consecutively to each of the equations has good smoothing properties for multigrid. The details of the relaxation have been described in Appendix B. The asymptotic smoothing rate, for the cases considered in this paper, is approximately 0.7 [35] which ensures that a multigrid coarse grid correction cycle with only a

few relaxations per gridlevel gives a grid independent convergence rate. This can be improved further by applying the relaxation in a line-wise manner. When the material is near incompressible ðν-0:5Þ the approach outlined in [43] needs to be used. C.2. Intergrid transfers and coarse grid operator The correction solved on the coarse grid can be interpolated to the fine grid using a standard linear interpolation, which can most efficiently be implemented one direction at the time. For the restriction of the residuals IH h the transpose of this interpolation can be used, referred to as full weighting. For the coarse grid operator the same discrete operator as on the fine grid can be used when the changes in the values of λ and μ are moderate. For cases with strong variations over small distances, a more sophisticated and robust approach is needed, see [35]. The coarsening is done for each variable separately. For each equation the coarse grid problem to be solved and the correction obtained from the coarser grid should be formulated as prescribed by the Full Approximation Scheme, i.e. Eqs. (35) and (37). This not only applies to the interior equations but also to the boundary conditions and to the force balance condition. This implies that each equation only has a zero right-hand side on the fine grid. On the coarser grids it has a righthand side that consists of the two parts given in (35). In the coarsening and refinement the boundary planes need to be treated separately from the interior. In these planes all transfers should be carried out in a 2D manner to and from the equivalent plane on the coarser grid. For the corrections a linear interpolation (in each dimension) for the corrections, and its transpose (full weighting) for the residuals, suits well. In the plane where contact occurs the transfers need to account for the complementarity of the contact problem see [15,19], by applying injection for the residuals of the normal stress equation when one of the points involved in the full weighting is a zero pressure point, and full weighting otherwise. This approach is identical to the approach for the elasto-hydrodynamic lubrication and non-adhesive contact problems used in [15,19]. C.3. Force balance This global constraint determines the mutual approach δ between the indenter and the stressed volume. The equation is relaxed according to (25). However, as changing the approach has a large scale effect on the entire solution, in the multigrid cycle it is either done in between cycles, or inside a cycle but then only on the coarsest grid. In the latter case it is essential that the coarsening of this condition is also formulated in the appropriate way, i.e. that the coarse grid formulation is constructed according to the Full Approximation Scheme equations, i.e. with a non-zero right-hand side which is obtained from the residual on the fine grid. Details of the implementation can be found in the example codes given in [19]. In either case, the change made to mutual approach slow down and even stall the convergence rate of the cycle when they are to big. When the changes are too small many cycles will be needed to achieve a solution for which the force balance condition is sufficiently satisfied, which deteriorates efficiency. C.4. Performance As an example of the performance of the developed method Fig. 12 shows the residual of the displacement equation in the z direction, and of the force balance equation as a function of the number of cycles for the most extreme case considered, the polycrystalline problem. Vð2; 1Þ cycles were used and the error

H. Boffy, C.H. Venner / Tribology International 92 (2015) 233–245

Fig. 12. Reduction of the residual of the z displacement equation and of the force balance equation as a function of the number of V ð2; 1Þ cycles for the polycrystalline case.

reduction is a factor 2 per cycle, which implies an effective error reduction per relaxation that is actually done of 0.8 which is quite close to the smoothing factor of the relaxation. Accounting for the complexity of the material topology and the effect of the force balance equation this performance is very satisfactory. Further improvement of the convergence rate can be achieved by applying line relaxation in each equation. Also a minor improvement can be achieved using W cycles. Concerning the computational cost, the calculations were done on a computer with a single processor. A CPU time of about 2.5 min is required to perform 1 V-cycle involving 30 million gridpoints (around 90 millions of unknowns for the three equations). On a finer grid with half the mesh size and about 230 million gridpoints (700 million unknowns) this time increases to about 20 min, consistent with the fact that O(N) complexity is predicted. Depending on the accuracy the user is looking for, the amount of cycles to perform may vary, see Fig. 12. However, using a FMG solver as shown in Fig. 1 typically two cycles per level will be needed to obtain a solution with an error that is well below the discretization error that is made anyway. References [1] Hughes TJR, Taylor RL, Sackman JCA, Kanoknukulchai W. A Finite Element Method for a class of contact impact problems. Comput Methods Appl Mech Eng 1976;8:249–76. [2] Hughes Jr T, Taylor RL, Sackman J. Finite element formulation and solution of contact impact problems in continuum. Technical report. Berkeley: University of California; 1978. [3] Curnier A. A theory of friction. Int J Solids Struct 1984;20:637–47. [4] Simo JC, Wriggers P, Taylor RL. A perturbed Lagrangian formulation for the Finite Element solution of contact problems. Comput Methods Appl Mech Eng 1985;50:163–80. [5] Wriggers P, Simo JC, Taylor RL. Penalty and augmented Lagrangian formulations for contact problems. In: Middleton J, Pande Balkema GN, editors. Proceedings of NUMETA, vol. 85. Swansea: British Technique Press; 1985. p. 97–105. [6] Bathe KJ, Chaudhary AB. A problems method for planar and axisymmetric contact problems. Int J Numer Methods Eng 1985;21:65–88. [7] Wriggers P. Computational contact mechanics. Berlin, Heidelberg: Springer; 2006. [8] Schrdöer A. Mixed Finite Element Methods of higher-order for model contact problems. SIAM J Numer Anal 2011;49:2323–39. [9] Hild P, Renard Y. An improved a priori error analysis for finite element approximations of Signorini's problem. SIAM J Numer Anal 2012;50:2400–19. [10] Liu T, Liu G, Wang QJ. An element-free Galerkin-Finite Element coupling Method for elasto-plastic contact problems. J Tribol 2006;128:1–9. [11] Oysu C, Fenner RT. Coupled FEM-BEM for elastoplastic contact problems using Lagrange multipliers. Appl Math Model 2006;30:231–47. [12] Hertz H. Über die berührung fester elastischer körper. J Reine Angew Math 1881;92:156–71.

245

[13] Kalker JJ, Allaert HJC, de Mul J. Nonlinear Finite Element Analysis in structural mechanics, Part V. In: The numerical calculation of contact problem in the theory of elasticity. Berlin, Heidelberg: Springer; 1981. p. 637–54. [14] Brandt A, Lubrecht AA. Multilevel matrix multiplication and fast solution of integral equations. J Comput Phys 1990;90:348–70. [15] Lubrecht AA, Ioannides A. A fast solution to the dry contact problem and the associated subsurface stress field, using multilevel techniques. J Tribol 1991;113:128–33. [16] Polonsky I, Keer L. A numerical method for solving rough contact problems based on the Multi-Level Multi-Summation and Conjugate Gradients techniques. Wear 1999;231:206–19. [17] Sainsot P, Lubrecht AA. Efficient solution of the dry contact for rough surfaces: a comparison of Fast Fourier Transform and Multigrid methods. Proc Inst Mech Eng Part J: J Eng Tribol 2011;225:441–8. [18] Vollebregt EAH. A new solver for the elastic normal contact problem using conjugate gradients, deflation, and an FFT-based preconditioner. J Comput Phys 2014;257:333–51. [19] Venner CH, Lubrecht AA. Multi-level methods in lubrication. Elsevier tribology series, Amsterdam; 2000. [20] Vorovich II, Ustinov IA. Pressure of a die on an elastic layer of finite thickness. J Appl Math Mech 1959;23:445–55. [21] Keer LM. The contact stress problem for an elastic sphere indenting an elastic layer. J Appl Mech 1964;31:143–5. [22] Kennedy FE, Ling FF. Elasto-plastic indentation of a layered medium. ASME J Eng Mater Technol 1974;96:97–103. [23] Plumet S, Dubourg MC. A 3-D model for a multilayered body loaded normally and tangentially against a rigid body: application to specific coatings. J Tribol 1998;120:668–76. [24] Polonsky I, Keer L. A fast and accurate method for numerical analysis of elastic layered contacts. J Tribol 2000;122:30–5. [25] Cai S, Bhushan B. A numerical three-dimensional contact model for rough, multilayered elastic/plastic solid surfaces. Wear 2005;259:1408–23. [26] Watremetz B, Baietto-Dubourg MC, Lubrecht AA. 2D thermo-mechanical contact simulations in a functionally graded material: a multigrid-based approach. Tribol Int 2007;40:754–62. [27] Boffy H, Baietto MC, Sainsot P, Lubrecht AA. An efficient 3D model of heterogeneous materials for elastic contact applications using Multigrid methods. J Tribol 2012;134:021401 DOI:101115/14006296. [28] Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc R Soc Lond 1957;A241:376–96. [29] Eshelby JD. The elastic field outside an elastic inclusion. Proc R Soc Lond 1959; A252:561–9. [30] Leroux J, Fulleringer B, Nelias D. Contact analysis in presence of spherical inhomogeneities within a half-space. Int J Solids Struct 2010;47:3034–49. [31] Koumi KE, Zhao LV, Leroux J, Chaise T, Nelias D. Contact analysis in the presence of an ellipsoidal inhomogeneity within a half space. Int J Solids Struct 2014;51:1390–402. [32] Zhou K, Chen WW, Keer LM, Ai X, Sawamiphakdi K, Glaws P, et al. Multiple 3D inhomogeneous inclusions in a half space under contact loading. Mech Mater 2011;43:444–57. [33] Zhou K, Hoh HJ, Wang X, Keer LM, Pang JHL, Song B, et al. A review of recent works on inclusions. Mech Mater 2013;60:144–58. [34] Zhou K, Wei R. Modeling cracks and inclusions near surfaces under contact loading. Int J Mech Sci 2014;83:163–71. [35] Boffy H, Venner CH. Multigrid solution of the 3D stress field in strongly heterogeneous materials. Tribol Int 2014;74:121–9. [36] Morales-Espejel GE, Boffy H, Venner CH, Ehret P. On the feasibility of the prediction of effects of material heterogeneity on surface fatigue for dry and lubricated contacts in realistic rolling bearing contacts. Tribol Int 2015, submitted for publication. [37] Briggs WL, Henson VE, McCormick SF. A multigrid tutorial. SIAM; Philadelphia, USA, 2000. [38] Trottenberg U, Oosterlee CW, Schueller A. Multigrid. London, UK: Academic Press; 2000. [39] Holmberg K, Matthews A. Coatings tribology, properties, mechanisms, techniques and applications in surface engineering. Tribology and interface engineering series, no. 56. Elsevier; Amsterdam, The Netherlands, 2009. [40] Stevanovic M, Yovanovich MM, Culham JR. Modeling contact between rigid sphere and elastic layer bonded to rigid substrate. IEEE Trans Compon Packag Technol 2001;24:207–12. [41] Chen WT, Engel PA. Impact and contact stress analysis in multilayer media. Int J Solids Struct 1972;8:1257–81. [42] Lele SK. Compact Finite Difference schemes with spectral-like resolution. J Comput Phys 1992;103:16–42. [43] Zhu Y, Sifakis E, Teran J, Brandt A. An efficient Multigrid Method for the simulation of high-resolution elastic solids. ACM Trans Graph 2010;29:1–18. [44] Brandt A. Multi-level adaptative solutions to boundary-value problems. Math Comput 1977;31:333–90. [45] Brandt A. Multigrid techniques: 1984 guide with application to fluid dynamics. Bonn: Gesellschaft für Mathematik und Datenverarbeitung; 1984.