- Email: [email protected]
Wear simulation using an incremental sliding Boundary Element Method
Wear 260 (2006) 1119–1128
Wear simulation using an incremental sliding Boundary Element Method G.K. Sfantos, M.H. Aliabadi ∗ Department of Aeronautics, Faculty of Engineering, Imperial College, University of London, South Kensington Campus, London SW7 2AZ, UK Received 20 January 2005; received in revised form 18 July 2005; accepted 28 July 2005 Available online 16 September 2005
Abstract In this paper, an application of the Boundary Element Method (BEM) to wear simulation problems is presented. The BEM is used for modelling the two bodies in contact. The material loss of the bodies is represented in terms of the applied load and the sliding distance and modelled using a linear wear model. An incremental method is adopted to solve the non-linear wear problem. The analysis demonstrates that the BEM is efficient for modelling the bodies in contact. The computational time remains low as only variables at the boundaries are calculated instead of the whole domain. The size of every sliding increment must be optimized in order the total CPU time to remain low and the final geometries of the bodies to be smooth, in some specified tolerances. © 2005 Elsevier B.V. All rights reserved. Keywords: Boundary Element Method; Wear modelling; Contact analysis
1. Introduction In most engineering applications involving multi-body contact wear can occur. The selection of appropriate materials and the suitable design of the contact bodies can be used for decreasing the wear. From the design point of view, the precise study and knowledge of the stress-strain fields developed in the components close to the contact region is required, as wear is closely related to the developed bearing pressure inside the contact area. The wear depth is related to the contact pressure distribution and as bodies lose material the distribution is altered. The pressure distribution is a function of the contact area and the contact area changes with respect to the sliding distance. To simulate these changes it is necessary to divide the total sliding distance into small increments and the whole problem to be solved for every sliding increment, updating the geometries of the bodies at every step. That is the Incremental Method which is usually used for standard wear simulation problems [1,2], for simulating machining processes [3,4] and for simulating wear in bioengineering as in total hip and knee prosthesis [5–7]. The solution converges as the sliding increment decreases.
∗
Corresponding author. Tel.: +44 20 7594 5077; fax: +44 20 7594 5078. E-mail address: [email protected] (M.H. Aliabadi).
0043-1648/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.wear.2005.07.020
One of the most popular numerical methods for modelling wear problems has been the Finite Element Method (FEM) [1–7], due to its ability to model complicated geometries that come into contact and thus wear out. An alternative numerical method is the Boundary Element Method (BEM). The BEM has enjoyed some popularity for modelling contact problems as from the numerical point of view the contact is on the boundary of the bodies, and therefore a boundary-type rather than a domaintype solution procedure, would be more suitable to the analysis of these types of problems. The application of the BEM to simple wear problems was reported by Serre et al. [8]. In their work, a conformal contact problem was modelled as a one-body-wear formulation. In this paper, a Boundary Element Method is presented for modelling wear in multi-body contact. The method presented here is capable for modelling conformal and non-conformal contact problems with one or more bodies losing material. Several examples are presented to demonstrate the accuracy and robustness of the proposed method.
2. Contact mechanics The first applications of the BEM to contact problems can be traced back to 80’s [9–11]. The formulation used in the present paper for modelling two solids in contact is similar to
1120
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
nodes of the potential contact-worn area of each body, which is ¯H ¯H defined to be the same for both bodies. The matrices H c and Gc H are of size 2M × 2Mc and their terms correspond to the local coordinates of each node of the potential contact-worn area. In H H addition, the uH nc and tnc have 2(M − Mc ) components. Both of them contain field unknowns and prescribed boundary values of the non-contact boundaries, in terms of displacements and trac¯H tions, respectively. The vectors u¯ H c and tc have 2Mc components and both contain the unknown boundary conditions, in terms of displacements and tractions, respectively, of the potential contact area, in local coordinates. To correlate the field unknowns inside the contact area of both bodies in contact, the contact boundary conditions must be introduced. These conditions relate the displacements and the tractions, in local coordinates, of both bodies. Therefore, to H transform the HH c and Gc matrices from global to local coordinates the following transformation matrix RH is used.
Fig. 1. The boundary division to contact and non-contact region.
⎡ the Multi-Region formulation [12]. In this section, a brief review of the boundary integral equation formulation for modelling the multi-body contact, is presented. For more details the authors are referred to [13,14]. The potential contact-wear boundary is denoted as ScH and H , for each body H = {A, B}, as the non-contact boundary as Snc Fig. 1 illustrates. H S H = Snc ∪ ScH .
(1)
The displacement integral equation [12] for the bodies shown in Fig. 1, can be written as:
nix
⎢ −ni ⎢ y ⎢ ⎢ H R = ⎢ ... ⎢ ⎢ ⎣ 0 0
niy
···
0
nix .. . 0
···
0 .. .
···
c nM x
0
···
c −nM y
..
.
0
⎤H
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ M c ny ⎦ 0 .. .
(4)
c nM x
for i = 1, Mc , where nx and ny are the components of the outward unit normal vector of each contact node, of each body H = {A, ¯H ¯H B}. The matrices H c and Gc are calculated using the following expressions:
(2) H where uH j , tj are components of displacements and tractions, respectively, for each body H = {A, B}, UijH (x , x), TijH (x , x) are the fundamental solution representing the displacements and the tractions, respectively, in the j direction at a point x due to a unit point force in the i direction at point x (Appendix A), and CijH (x ) is the so-called free term. H of each body To solve Eq. (2) the boundaries ScH and Snc H H H = {A, B} are discretized into Nc and Nnc elements, respectively. Each potential contact element is composed by mH c number of nodes and each non-contact element by mH . After the nc discretization and using the point collocation method for solution, Eq. (2) can be written in matrix form as uH tH nc nc H H H H ¯ ¯c = Gnc G (3) Hnc H H c ¯ u¯ H t c c
for each body H = {A, B}. H H H The matrices HH nc and Gnc are of size 2M × 2(M − Mc ), for two dimensional bodies, containing known integrals of the product of the shape functions, the Jacobians and the fundamental fields TijH and UijH , respectively. MH is the total number of nodes composing each body H = {A, B} and Mc is the total number of
H H ¯H H c = R HC ¯ H = RH GH . G c
(5)
C
In the case of a frictionless contact, where the coefficient of friction is either too low or practically zero, the contact modes for a non-conformal contact are reduced to Slip and Separation only [15,13]. The same situation occurs in the case of a sliding contact, where all the contact node pairs slide over each other. Additionally as friction does not affect much the normal pressure distribution [14], the coefficient of friction can be neglected for simplicity, as in the frictionless case. Therefore, considering a sliding wear non-conformal problem, where one body slides over another, the contact modes of the contact node pairs are Slip or Separation only. The boundary conditions implemented for the aforementioned contact modes and the check criteria for any contact mode violations are presented in Appendix B. According to the Multi-Region technique [12], the individual systems of equations are combined for the two (or more) regions and the interface boundary conditions are subsequently applied [13,14]. Considering Eq. (3) for each body H = {A, B}, and after
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
rearranging and applying the pre-scribed boundary conditions of the non-contact boundary, the resulting system of equations can be written as:
tion Wear, which arises from the adhesive forces set up whenever atoms come into intimate contact [17,18]. In delamination wear, the sliding speed is low enough that surface heating can be neglected. The model that describes this mechanism is the well known Archard’s linear wear equation [19], firstly proposed by Holm [20], given as: W H = KH
(6) H where the matrices AH nc and Rnc contain terms of both the H H previous matrices Hnc and Gnc , after the substitution of the prescribed boundary conditions of the non-contact area of each H body H = {A, B}. The column vectors xH nc and ync contain the unknown and the known, respectively, components of tractions and displacements of the non-contact boundary, after the substitution of the pre-scribed boundary values. The submatrix BC1 in Eq. (6) has 2Mc × 8Mc components and contains all the boundary conditions for the displacements and the tangential tractions for the contact node pairs in slip mode and also the constraints for the zero tangential and normal tractions for the pairs in separation mode. The submatrix BC2 has 2Mc × 4Mc components and contains all the boundary conditions for the tractions equilibrium inside the contact-worn area. The column matrices F1 and F2 have 2Mc components each and contain the right-hand sides of the boundary conditions implemented in sub-matrices BC1 and BC2 (Appendix C). The above coupled system, containing the elastostatic problems of both bodies plus the contact problem, is square of size Mt × Mt , where Mt = 2(MA + MB ) + 4Mc . Taking a closer look at Eq. (6) it can be seen that the square matrix is highly sparsed. More than 60% of the total components are zero terms. This becomes significant as wear simulation demands solution of the above system several times, after updating the geometry of the worn body at every step of the solution algorithm. Another advantage of the specific method is that until convergence to be achieved concerning the modes of contact for each node pair, only some components of the submatrix BC1 on the left-hand side of Eq. (6) may change. Moreover as the changes take place only in a small part of the overall system, there is no need to calculate every time the inverse of the total matrix; however, using solution algorithms such as Sherman–Morrison or Woodbury formulas [16] the final solution of the above system can be speeded up several times [12,13].
3. Wear model In the case where the sliding speed remains at low levels and the applied load does not exceed a limit where seizure takes place, the wear that occur is characterized as Delamina-
1121
P S Ho
(7)
where WH is the volumetric material loss of the body in consideration, H = {A, B}, KH the wear coefficient, Ho the hardness of the softer body in contact, P the applied load and SH the sliding distance the body has travelled. For any pair of materials that slide over each other, a specific constant KH can be determined through an experiment. Using a wear model in a computational technique for wear simulation, the parameter “sliding material” is directly introduced in the formulation. However, the model described by Eq. (7) does not provide sufficient details; for that reason instead of the wear volume the Wear Depth can be used for modelling the worn geometry. Moreover, the wear depth is of practical importance in every day engineering as it can be measured directly on a specific mechanical component. In the two-dimensional case, the above model describes the volume of material loss per unit thickness, with respect to the sliding distance, when a specific load P per unit thickness is applied. Therefore, differentiating Eq. (7) with respect to the width of the bearing surface, results to: hH (x, s) = kH p(x, s)sH
(8)
where hH is the wear depth on point x of body H = {A, B} after sH sliding distance, p the pressure distribution developed inside the contact area and kH the modified wear coefficient given by kH = KH /Ho . Both pressure and wear depth are functions of the location and the sliding distance as the bodies lose material during the wear process and the topography of their bearing surfaces changes. Furthermore, in cases as a pin on a rotating disc problem, the travelled distance sH by each point x on the bearing surfaces is also a function of the location of the point. However, this sliding distribution with respect to the location of each point on the bearing surfaces can be evaluated directly by the kinematics of the specific tribosystem. The direction of wear, the wear vector, is taken to coincide with the pressure vector; that is normal to the bearing surface of the body in contact. 4. Solution algorithm The most commonly used technique for simulating wear problems is the Incremental Method [1–7]. It is clear from Eq. (8) that once the pressure distribution inside a contact area is known for a specific sliding distance, the wear depth, with respect to the width of the contact area, can be calculated. However, as the function of the pressure distribution is not known a priori (due to the variation with the sliding distance), the wear depth cannot be calculated directly. To overcome this difficulty the
1122
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
geometry, where the Boundary Element mesh (or the FE mesh) is generated. Then the contact analysis starts with an initial guess of the contact node pairs. The coupled contact problem is solved and checks for contact mode violations are performed. If the initial guess for the contact node pairs that come into contact, for the specific applied normal load, was false, the assumed pairs are updated and the problem is solved again and again until no contact mode violations occur. Then using the resulting pressure distribution inside the contact area, and using Eq. (8), the resulting wear depth is calculated. The geometries of both bodies are updated (in case that both bodies lose material), and the intermediate results are stored. If the accumulated sliding distance reaches the total sliding distance then the simulation is finished. Otherwise a new BE mesh is generated and the above procedure is repeated. Some fine tuning of this method is required in terms of the sliding-time increment that will be used. If the increment is relatively large, then oscillations will occur to the resulting geometries and consequently will not be smooth. Moreover, the validity of the calculated wear depth will not be definite. A trial and error procedure is adopted to calibrate the maximum step that produces smooth geometries and wear depth, and additionally the computational time to be reasonable small. Another way is to use a variable step technique, in order the maximum wear depth for every increment to be less than a specified value [1]; thus, the solution procedure may speed up. 5. Results and discussion 5.1. One body wear formulation
Fig. 2. Solution algorithm for the Incremental Method.
total sliding distance is incremented and the problem is solved for small increments updating every time the geometries of the bodies. The basic algorithm of the Incremental Method using the BEM is illustrated in Fig. 2. The algorithm begins with the initial
In the case where only one of the two bodies in contact loses material, the procedure to simulate a wear experiment is relatively straightforward. To validate the BEM presented here, analytical results using the wear model and the results from the work of P˜odra and Andersson [1] were used. P˜odra and Andersson [1] used the Finite Element Method ANSYS code. In their work 3D formulation was used for simulating the classic problem of a pin on disc wear experiment. Their results were compared, in terms of the wear depth, with the experimental evidences that were used to calibrate the wear model. In the present work, a two-dimensional formulation is employed. For this reason and to allow comparison with the results in [1], the calculated wear depth was transformed to wear volume, using the geometric equation that calculates the volume for a specific depth of a spherical ended pin, given by [21]: WA =
π A 2 A 3R − hA hm m . 3
(9)
For the maximum wear depth hA m calculated using the FEM package, and observed experimentally in the work of P˜odra and Andersson [1], and for a cylindrical ended pin of thickness 1 mm,
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128 Table 1 Parameters used in example number 1 Wear coefficient (Pa−1 ) Total sliding distance (m) Normal load (N/mm) Poisson ratio Radius of the pin (m) Elastic moduli (GPa) Sliding increment (mm)
kA = 1.33 ± 0.54 × 10−13 S = 3.0 P = 50 vA = vB = 0.3 RA = 5 × 10−3 EA = EB = 210 ds = 5, 10, 20, 40, 50
given by [21]: RA − hA m W = (R ) cos RA A 2 − RA − hA 2RA hA m − (hm ) m A
A 2
−1
(10)
for the maximum wear depth hA m calculated using the BEM in the present work. The reason that wear volume is used instead of wear depth to compare the two different formulations (i.e. a 2D and a 3D), is due to the wear law, Eq. (7), as it is not taking into account the geometries of the bodies. Hence, using the same normal applied load, wear coefficient and sliding distance, the wear volume must be exactly the same whatever are the contacting geometries. Each body was modelled using 65 quadratic elements, whereas 69 nodes (34 elements of the 65) were used to model the potential contact-worn area. The parameters of the specific problem are presented in Table 1. The mesh density remains constant throughout the process and consequently the first steps are more crucial since they involve less number of nodes to simulate high
1123
peak contact pressures. As wear progresses, smoother pressures will result and more nodes of the potential contact-worn surfaces will come into contact. Fig. 3 illustrates the initial pressure distribution, before wear, evaluated analytically by Hertzian theory [15] and numerically by the BEM (plane strain). Several trials considering different mesh sizes were carried out to assure the accurate evaluation of the bearing pressure. The mesh size used in the present example (69 nodes for the potential contact-worn area) corresponds to part (c) of Fig. 3, incorporating 9 nodes inside the initial contact area. It is obvious that the developed bearing pressure is evaluated sufficiently accurate, even during the initial step where only 9 nodes are in contact. The same parameters were used in the work of P˜odra and Andersson [1], who calibrate the wear coefficient kA experimentally. In this problem, only the pin is assumed to lose material, as the disc is assumed very hard and its wear can be neglected compared to that of the pin. Moreover, the disc diameter is assumed very large compared to the pin diameter and the position of the pin is far from the center of revolution of the disc. Consequently, as the disc rotates, the difference of the travelled distance between points on the bearing surface of the pin can be neglected and thus all points on the pin, either in contact or not, is assumed to travel exactly the same distance [1]. In case where the travelled distance between points corresponding to the internal and the external radius of the sliding path varies significantly, then the sliding distance distribution must be evaluated from the kinematics of the specific tribosystem, corresponding to every node of the discretized geometry of the body that loses material. The modified wear law, Eq. (8), is applied to all points on the bearing surface of the pin; for points that are in contact and consequently under bearing pressure, the corresponding wear
Fig. 3. Initial contact pressure: (a) #3 nodes in contact; (b) #5 nodes in contact; (c) #9 nodes in contact; (d) #13 nodes in contact (a: width of contact area).
1124
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
Fig. 4. Comparison of the wear volume calculated analytically by the wear model with the wear volume evaluated using the BEM and the Incremental Method, with respect to the sliding distance.
depth is calculated by Eq. (8); however, for points that are not in contact yet, the corresponding wear depth is zero as these points are not under bearing pressure. Therefore, as the pin slides over the disc, the travelled distance by each point on the pin’s bearing surface coincide with the sliding distance of the pin (as a body), while each point on the pin’s bearing surface slides different distance as it comes into contact at different time as the pin loses material. In this paper, the following analysis corresponds to the sliding distance of the pin, as a body. Fig. 4 illustrates the calculated wear volume by the present work, using Eq. (10) in comparison with the analytically calculated wear volume using directly the wear model, Eq. (7) after substituting kH = KH /Ho , and by using the same parameters presented in Table 1. It is apparent that as the sliding increment decreases the wear volume converges. It is obvious that the numerically evaluated wear volume, using the BEM and the incremental form of the wear model, are in good agreement compared with the analytically calculated wear volume. Additionally, in Fig. 5, the maximum wear depth calculated by the proposed formulation is illustrated against the wear depth calculated solving inversely Eq. (10) and considering the wear volume that was evaluated using the wear model, Eq. (7), which
Fig. 6. Comparison of the wear volume evaluated experimentally, by FEM and BEM, with respect to the sliding distance.
is illustrated in Fig. 4. The same behaviour is demonstrated as the sliding increment decreases. Fig. 6 illustrates the wear volume calculated in the present work, compared to the wear volume from [1], by using Eq. (9) for both the numerical and the experimental results that were reported in [1]. The upper and lower limits illustrated in the same figure are the deviation of the wear coefficient k as it was calibrated experimentally by P˜odra and Andersson [1]. In Table 2, the parameters of another numerical example are presented. In this example each body was modelled using 50 quadratic elements in total, whereas 49 nodes (24 elements) were used to compose the potential contact area. Fig. 7 illustrates the final geometry of the worn pin after a total sliding distance of 500 mm, while Fig. 8 illustrates the deviation of the maximum wear depth for various sliding increments. It is obvious from both figures that as the sliding increment decreases the solution converges. However the computational time increases almost linearly as the same calculation processes must be carried out for a greater number of sliding increments. Taking a closer look to Fig. 7, it can be seen that the use of large sliding increments results to the development of an artificial roughness on the bearing surface of the body that wears out. This is due to the incrementation of the sliding distance, as large increments create large craters which produce localized high pressures at their edges after the geometry update at the next step. These localized high pressures create new large craters and the phenomenon continues. As a result geometries with artificial roughness are created. On the other hand, as the sliding increment decreases, these craters are relatively small Table 2 Parameters used in example number 2
Fig. 5. Variation of analytically and numerically calculated maximum wear depth with respect to the sliding distance.
Wear coefficient (Pa−1 ) Total sliding distance (m) Normal load (N/mm) Poisson ratio Radius of the pin (m) Elastic moduli (GPa) Sliding increment (mm)
kA = 1.33 × 10−13 S = 0.5 P = 4.63 vA = vB = 0.3 RA = 10 × 10−3 EA = EB = 70 ds = 2.5, 5, 10, 20
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
1125
Fig. 7. Deviation of the final geometry with respect to the sliding increment: (a) ds = 20 mm; (b) ds = 10 mm; (c) ds = 5 mm; (d) ds = 2.5 mm (a: width of contact-worn area).
and produce a smooth transition from the unworn to the worn geometry. The aforementioned artificially created roughness is also responsible for the deviation of the wear depth and volume, as Figs. 4, 5 and 8 illustrate. 5.2. Two bodies wear formulation In cases where both bodies wear out in comparable time, the problem becomes more complicated as both bearing surfaces change with respect to the sliding distance. In such cases, it may be of interest to simulate the evolution of wear on each body in contact individually. To do so the modified wear law, Eq.
(8), must be used for each body H = {A, B}. Table 3 illustrates the parameters used in the case of two bodies wear formulation. In the specific problem it is assumed that all points on the pin’s bearing surface travel the same distance with all the points on the disc’s bearing surface. This is done for simplicity and to demonstrate the relation between the final wear of each body and its wear coefficient. In practise, in a pin on a rotating disc problem, the sliding distance distribution must be evaluated, corresponding to the changeable location of every node on both bearing surfaces, considering the kinematics of the tribosystem and the contact duration, as it is already explained in Section 5.1. Fig. 9 illustrates the final geometry of the two bodies. Taking a closer look to the specific figure, it can be seen that the maximum wear depth of the pin, body A, is exactly two times the maximum wear depth of the disc, body B. The specific relation between the maximum wear depth of the two bodies results from the input data, Table 3, where the wear coefficient of the pin is exactly two times the coefficient of the disc. In the case where both bodies had the same wear coefTable 3 Parameters for the two-bodies-wear formulation example
Fig. 8. Deviation of the maximum wear depth with respect to the sliding increment.
Wear coefficient (Pa−1 ) Total sliding distance (m) Normal load (N/mm) Poisson ratio Radius of the pin (m) Elastic moduli (GPa) Sliding increment (mm)
kA = 0.887 × 10−13 , kB = 0.443 × 10−13 S = 0.5 P = 4.63 vA = vB = 0.3 RA = 10 × 10−3 EA = EB = 70 ds = 2.5
1126
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
In Fig. 11, the overall maximum wear depth of a wear problem where both bodies lose material, is compared with a problem where only one body lose material, but its wear coefficient is equal to the sum of the two coefficient of the two bodies wear formulation problem. It is obvious that the resulting overall wear depth from the two bodies wear formulation is the same with the one body wear formulation, as the summation of the wear coefficients in the first case is equal to the wear coefficient in the second case. This analysis demonstrates the linear behaviour of the incremented modified wear law, Eq. (8), with respect to the wear coefficients kH . 6. Conclusions Fig. 9. Two bodies wear formulation: Final geometries for k(pin) = 2k(disc) (a: width of contact-worn area).
Fig. 10. Two bodies wear formulation: final geometries for k(pin) = k(disc) (a: width of contact-worn area).
ficients, i.e. kA = kB = 0.665 × 10−13 Pa−1 , the resulting geometries are illustrated in Fig. 10. In this figure the maximum wear depth of the pin is equal to the maximum wear depth of the disc.
In the present work, a two-dimensional application of the Boundary Element Method in wear simulation was presented. The Incremental Method was used where the total sliding distance is divided into small increments. The concluding remarks of the present work are summarized as follows: (1) The stress-strain fields inside the contact area are evaluated fast and accurately using the BEM. Thus, monitoring of the changes in these fields during or after wear, results to a better knowledge of the wear evolution on a specific component. (2) The computational time remains low as only boundary modelling is required in contrast to the FEM. This is a great advantage considering the wear simulation problems, where the geometry is updated after every step. Thus remeshing of the bodies in contact becomes relatively easy and fast. (3) The wear depth is calculated more accurately using the BEM, in comparison with the FEM, as it is directly related to the pressure distribution inside the contact. This is due to the ability of the BEM to evaluate with high accuracy both stresses and displacements unknowns on the boundaries of the bodies in concern. (4) The size of every sliding increment must be optimized in order the final geometries of the bodies to be smooth, in some specified tolerance, and the computational time to remain low. Appendix A The fundamental solutions in Eq. (2) can be obtained by several methods. Here, the fundamental solutions obtained using the Galerkin vector are used [12]. The displacement and traction fundamental solutions for the two-dimensional plain strain condition are given as follows:
1 1 Uij (x , x) = (3 − 4v) ln δij + r,i r,j 8πG(1 − v) r Tij (x , x) =
Fig. 11. Comparison of the maximum wear depth between one and two bodies wear formulation.
∂r [(1 − 2v)δij + 2r,i r,j ] ∂n −(1 − 2v)(r,i nj − r,j ni )
−1 4π(1 − v)r
where Uij (x , x) and Tij (x , x) represents the displacement and the traction, respectively, in the j direction at point x due to a unit
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
point force acting in the i direction at x . G is the shear modulus, nj denotes the components of the outward normal unit vector at the field point x, r is the distance between the field and the force point and δij is the Kronecker delta function. Appendix B The boundary conditions considering separation and slip modes are given in Table 4, where (tcH )t and (tcH )n are components of tangential and normal tractions, respectively, and (uH c )n are components of normal displacements, of each body H = {A, B}. The normal gap between the potential contact node pairs is defined as gapA–B . As shown in Fig. 12 the initial normal gap gapA–B is defined as the normal distance between a node pair, and is calculated as follows [13]: gapA–B = |dy|¯nA nA y + |dx|¯ x
1127
Table 5 Contact mode check criteria Decision . . . assumption
Separation
Contact
Separation Contact
B A−B (uA c )n + (uc )n < gap (tcA or B )n ≥ 0
B A–B (uA c )n + (uc )n ≥ gap (tcA or B )n < 0
Calculating the normal gap in this way, the surfaces of the two bodies will come into contact after the specific gap is closed, either if the nodes are in touch or slide over each other. In order to obtain the correct contact area without any geometrical incompatibility and contact mode violations, the displacements and tractions for every node in the contact zone, must be checked against the criteria presented in Table 5.
(11) Appendix C
¯A for all potential contact node pairs, where n¯ A y and n x are the H components of the average normal unit vector n¯ , of each body H = {A, B}, which is calculated from the outward normal unit vectors of the specific pair of nodes nA , nB by n¯ A =
nA − nB = −n¯ B . |nA − nB |
(12)
Table 4 Boundary constraints Separation
Slip
(tcA )t − (tcB )t = 0 (tcA )n − (tcB )n = 0 (tcA )t = 0 (tcA )n = 0
(tcA )t − (tcB )t = 0 (tcA )n − (tcB )n = 0 (tcA )t = 0 B A–B (uA c )n + (uc )n = gap
Fig. 12. Definition of the average normal unit vector and the initial normal gap.
The sub-system in Eq. (6), corresponding to the sub-matrix BC2 , includes all the boundary constraints for the tractions equilibrium inside the contact area. It can be written as: A ¯t [BC2 ] cB ¯tc ⎤ ⎡ 1 0 0 0 . . . −1 0 0 0 ... ⎥ ⎢ 0 ...⎥ ⎢ 0 1 0 0 . . . 0 −1 0 ⎥ ⎢ ⎢ 0 −1 0 . . . ⎥ = {F2 } ⇒ ⎢ 0 0 1 0 . . . 0 ⎥ ⎥ ⎢ 0 0 −1 . . . ⎥ ⎢0 0 0 1 ... 0 ⎦ ⎣ .. .. .. .. .. .. .. .. . . .. . . . . . . . . . . ⎧ ⎫ ⎧ ⎫ (tcA )n,1 ⎪ ⎪ ⎪ ⎪ ⎪0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (tcA )t,1 ⎪ ⎪ ⎪ ⎪ ⎪ 0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ A ⎪ ⎪ ⎪ ⎪ (t ) ⎪ ⎪ ⎪ ⎪ 0 c n,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ A ⎪ ⎪ ⎪ ⎪ (t ) ⎪ ⎪ ⎪ 0 ⎪ ⎪ c t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ .. ⎪ ⎪ .. ⎪ ⎪ ⎪ ⎨ ⎨ ⎬ .⎬ . × = B ⎪ ⎪ ⎪ (tc )n,1 ⎪ ⎪0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ B ⎪ ⎪ (tc ) ⎪ ⎪0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ t,1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ B ⎪ ⎪ ⎪ ⎪ 0 ⎪ ⎪ ⎪ ⎪ (t ) c ⎪ ⎪ ⎪ n,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ B 0 ⎪ ⎪ ⎪ (tc )t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ . ⎪ ⎩ ⎪ ⎪ .. ⎭ ⎩ .. ⎪ ⎭ . where (tcA )n,1 is the component of the normal traction acting on node 1 of the contact area of body A. This component is related to the normal traction acting on node 1 of the contact area of body B through the boundary condition for the tractions equilibrium (Table 4). Additionally, (tcA )t,1 is the component of the tangential traction acting on node 1 of the contact area of body A and is related to the corresponding tangential traction (tcB )t,1 acting on the coupled node of body B.
1128
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
If the contact node pair 1 is in Slip and the contact node pair 2 is in Separation mode, the sub-system in Eq. (6) corresponding to the sub-matrix BC1 can be written as: ⎧ A⎫ u¯ c ⎪ ⎪ ⎪ ⎪ B⎪ ⎪ A A ⎨ ⎬ [BC1 ]
u¯ c = {F1 } ⇒ [BC1,1 ] A ⎪ ⎪ ¯tc ⎪ ⎪
⎪ ⎩ ¯B ⎪ ⎭
u¯ c u¯ Bc
tc
⎡
1
0
0
0
...
.. .
.. .
.. .
.. .
..
1
0
0
0
.. .
.. .
.. .
.. .
⎢0 0 0 0 ... 0 0 0 0 ⎢ ⎢0 0 0 0 ... 0 0 0 0 ⇒⎢ ⎢ ⎢0 0 0 0 ... 0 0 0 0 ⎣ .
+ [BC1,2 ]
⎧ A ⎫ (uc )n,1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (uA ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ c ⎪ ⎪ A t,1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (u ) ⎪ c n,2 ⎪ ⎪ ⎤⎪ ⎪ ⎪ A ... ⎪ ⎪ ⎪ (uc )t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ...⎥ ⎪ ⎪ . ⎪ ⎥⎪ ⎨ ⎬ . ⎥ . ...⎥ ⎥ (uB ) ⎪ ⎪ c n,1 ⎪ ...⎥⎪ ⎪ ⎪ ⎪ ⎦⎪ ⎪ (uBc )t,1 ⎪ ⎪ ⎪ .. ⎪ ⎪ ⎪ ⎪ . ⎪ ⎪ B ⎪ ⎪ (u ) ⎪ c n,2 ⎪ ⎪ ⎪ ⎪ ⎪ (uB ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ c t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ . ⎩ . ⎪ ⎭ .
⎡
0
0
0
0
... 0
0
0
0
.. .
.. .
.. .
.. .
..
.. .
.. .
.. .
.. .
⎢0 1 0 0 ... 0 0 0 0 ⎢ ⎢0 0 1 0 ... 0 0 0 0 +⎢ ⎢ ⎢0 0 0 1 ... 0 0 0 0 ⎣ .
⎧ A ⎫ (t ) ⎪ ⎪ ⎪ ⎪ cA n,1 ⎪ ⎪ ⎪ ⎪ ⎪ (tc )t,1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ A ⎪ ⎪ (t ) c n,2 ⎪ ⎪ ⎪ ⎤⎪ ⎪ ⎪ A ... ⎪ ⎪ ⎪ (tc )t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ...⎥ ⎪ . ⎪ ⎪ ⎥⎪ ⎨ ⎬ . ⎥ . ...⎥ ⎥ (t B ) ⎪ ⎪ c n,1 ⎪ ...⎥⎪ ⎪ B ⎪ ⎪ ⎦⎪ ⎪ (t ) ⎪ ⎪ ⎪ .. ⎪ ⎪ Bc t,1 ⎪ ⎪ . ⎪ ⎪ ⎪ (tc )n,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ B ⎪ ⎪ (t ) ⎪ ⎪ c t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ .. ⎪ ⎭ .
=
⎧ ⎫ gapA−B ⎪ ⎪ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 0 ⎪ ⎪ ⎪ ⎪ ⎨ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩
0 0 .. .
⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
.
¯tc ¯tBc
= {F1 }
References [1] P. P˜odra, S. Andersson, Simulating sliding wear with finite element method, Tribol. Int. 32 (1999) 71–81. [2] P. P˜odra, S. Andersson, Finite element analysis wear simulation of a conical spinning contact considering surface topography, Wear 224 (1999) 13–21. [3] T. MacGinley, J. Monaghan, Modelling the orthogonal machining process using coated cemented carbide cutting tools, J. Mater. Process. Technol. 118 (2001) 293–300. [4] Y.-C. Yen, J. Sohner, B. Lilly, T. Altan, Estimation of tool wear in orthogonal cutting using the finite element analysis, J. Mater. Process. Technol. 146 (2004) 82–91. [5] T.A. Maxian, T.D. Brown, D.R. Pedersen, J.J. Callaghan, A slidingdistance-coupled finite element formulation for polyethylene wear in total hip arthroplasty, J. Biomech. 29 (1996) 687–692. [6] T.A. Maxian, T.D. Brown, D.R. Pedersen, J.J. Callaghan, 3-Dimensional sliding/contact computational simulation of total hip wear, Clin. Orthop. Rel. Res. 333 (1996) 41–50. [7] J.S.-S. Wu, J.-P. Hung, C.-S. Shu, J.-H. Chen, The computer simulation of wear behaviour appearing in total hip prosthesis, Comput. Methods Prog. Biomed. 70 (2003) 81–91. [8] I. Serre, M. Bonnet, R.-M. Pradeilles-Duval, Modelling an abrasive wear experiment by the boundary element method, C.R. Acad. Sci. Paris 329 (2001) 803–808, S´erie II b. [9] Anderson, T., Boundary Elements in two-dimensional Contact and Friction, Diss No 85, Linkoping Institute of Technology, 1982. [10] F. Paris, J.A. Garrido, On the use of discontinuous elements in 2D contact problems, in: Proc. of Boundary Elements VII, Computational Mechanics Publications, Southampton, 1985, 13.27–13.39. [11] G. Karami, R.T. Fenner, A two-dimensional BEM method for thermoelastic body forces contact problems, in: Proc. of Boundary Elements IX, Computational Mechanics Publications, Southampton, 1987, 417–419. [12] M.H. Aliabadi, The Boundary Element Method Applications in Solids and Structures, 2, Wiley, London, 2002. [13] K. Man, M.H. Aliabadi, D.P. Rooke, BEM Frictional Contact Analysis: modelling considerations, Eng. Anal. Bound. Elem. 11 (1993) 77–85. [14] K. Man, M.H. Aliabadi, D.P. Rooke, BEM Frictional Contact Analysis: An Incremental Loading Technique, Comput. Struct. 47 (1993) 93–905. [15] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 2001. [16] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes in Fortran, in: The Art of Scientific Computing, second ed., Cambridge University Press, Cambridge, 1994. [17] S.C. Lim, M.F. Ashby, Wear-mechanism maps, Acta Metall. 35 (1) (1987) 1–24. [18] E. Rabinowicz, Friction and Wear of Materials, second ed., Wiley, New York, 1995. [19] J.F. Archard, Contact and rubbing of flat Surfaces, J. Appl. Phys. 24 (1953) 981–988. [20] R. Holm, Electric Contacts, Gebers, Stockholm, 1946. [21] M.R. Spiegel, J. Liu, Mathematical Handbook of Formulas and Tables, ISE ed., Mc-Graw Hill, New York, 1999.
Wear simulation using an incremental sliding Boundary Element Method G.K. Sfantos, M.H. Aliabadi ∗ Department of Aeronautics, Faculty of Engineering, Imperial College, University of London, South Kensington Campus, London SW7 2AZ, UK Received 20 January 2005; received in revised form 18 July 2005; accepted 28 July 2005 Available online 16 September 2005
Abstract In this paper, an application of the Boundary Element Method (BEM) to wear simulation problems is presented. The BEM is used for modelling the two bodies in contact. The material loss of the bodies is represented in terms of the applied load and the sliding distance and modelled using a linear wear model. An incremental method is adopted to solve the non-linear wear problem. The analysis demonstrates that the BEM is efficient for modelling the bodies in contact. The computational time remains low as only variables at the boundaries are calculated instead of the whole domain. The size of every sliding increment must be optimized in order the total CPU time to remain low and the final geometries of the bodies to be smooth, in some specified tolerances. © 2005 Elsevier B.V. All rights reserved. Keywords: Boundary Element Method; Wear modelling; Contact analysis
1. Introduction In most engineering applications involving multi-body contact wear can occur. The selection of appropriate materials and the suitable design of the contact bodies can be used for decreasing the wear. From the design point of view, the precise study and knowledge of the stress-strain fields developed in the components close to the contact region is required, as wear is closely related to the developed bearing pressure inside the contact area. The wear depth is related to the contact pressure distribution and as bodies lose material the distribution is altered. The pressure distribution is a function of the contact area and the contact area changes with respect to the sliding distance. To simulate these changes it is necessary to divide the total sliding distance into small increments and the whole problem to be solved for every sliding increment, updating the geometries of the bodies at every step. That is the Incremental Method which is usually used for standard wear simulation problems [1,2], for simulating machining processes [3,4] and for simulating wear in bioengineering as in total hip and knee prosthesis [5–7]. The solution converges as the sliding increment decreases.
∗
Corresponding author. Tel.: +44 20 7594 5077; fax: +44 20 7594 5078. E-mail address: [email protected] (M.H. Aliabadi).
0043-1648/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.wear.2005.07.020
One of the most popular numerical methods for modelling wear problems has been the Finite Element Method (FEM) [1–7], due to its ability to model complicated geometries that come into contact and thus wear out. An alternative numerical method is the Boundary Element Method (BEM). The BEM has enjoyed some popularity for modelling contact problems as from the numerical point of view the contact is on the boundary of the bodies, and therefore a boundary-type rather than a domaintype solution procedure, would be more suitable to the analysis of these types of problems. The application of the BEM to simple wear problems was reported by Serre et al. [8]. In their work, a conformal contact problem was modelled as a one-body-wear formulation. In this paper, a Boundary Element Method is presented for modelling wear in multi-body contact. The method presented here is capable for modelling conformal and non-conformal contact problems with one or more bodies losing material. Several examples are presented to demonstrate the accuracy and robustness of the proposed method.
2. Contact mechanics The first applications of the BEM to contact problems can be traced back to 80’s [9–11]. The formulation used in the present paper for modelling two solids in contact is similar to
1120
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
nodes of the potential contact-worn area of each body, which is ¯H ¯H defined to be the same for both bodies. The matrices H c and Gc H are of size 2M × 2Mc and their terms correspond to the local coordinates of each node of the potential contact-worn area. In H H addition, the uH nc and tnc have 2(M − Mc ) components. Both of them contain field unknowns and prescribed boundary values of the non-contact boundaries, in terms of displacements and trac¯H tions, respectively. The vectors u¯ H c and tc have 2Mc components and both contain the unknown boundary conditions, in terms of displacements and tractions, respectively, of the potential contact area, in local coordinates. To correlate the field unknowns inside the contact area of both bodies in contact, the contact boundary conditions must be introduced. These conditions relate the displacements and the tractions, in local coordinates, of both bodies. Therefore, to H transform the HH c and Gc matrices from global to local coordinates the following transformation matrix RH is used.
Fig. 1. The boundary division to contact and non-contact region.
⎡ the Multi-Region formulation [12]. In this section, a brief review of the boundary integral equation formulation for modelling the multi-body contact, is presented. For more details the authors are referred to [13,14]. The potential contact-wear boundary is denoted as ScH and H , for each body H = {A, B}, as the non-contact boundary as Snc Fig. 1 illustrates. H S H = Snc ∪ ScH .
(1)
The displacement integral equation [12] for the bodies shown in Fig. 1, can be written as:
nix
⎢ −ni ⎢ y ⎢ ⎢ H R = ⎢ ... ⎢ ⎢ ⎣ 0 0
niy
···
0
nix .. . 0
···
0 .. .
···
c nM x
0
···
c −nM y
..
.
0
⎤H
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ M c ny ⎦ 0 .. .
(4)
c nM x
for i = 1, Mc , where nx and ny are the components of the outward unit normal vector of each contact node, of each body H = {A, ¯H ¯H B}. The matrices H c and Gc are calculated using the following expressions:
(2) H where uH j , tj are components of displacements and tractions, respectively, for each body H = {A, B}, UijH (x , x), TijH (x , x) are the fundamental solution representing the displacements and the tractions, respectively, in the j direction at a point x due to a unit point force in the i direction at point x (Appendix A), and CijH (x ) is the so-called free term. H of each body To solve Eq. (2) the boundaries ScH and Snc H H H = {A, B} are discretized into Nc and Nnc elements, respectively. Each potential contact element is composed by mH c number of nodes and each non-contact element by mH . After the nc discretization and using the point collocation method for solution, Eq. (2) can be written in matrix form as uH tH nc nc H H H H ¯ ¯c = Gnc G (3) Hnc H H c ¯ u¯ H t c c
for each body H = {A, B}. H H H The matrices HH nc and Gnc are of size 2M × 2(M − Mc ), for two dimensional bodies, containing known integrals of the product of the shape functions, the Jacobians and the fundamental fields TijH and UijH , respectively. MH is the total number of nodes composing each body H = {A, B} and Mc is the total number of
H H ¯H H c = R HC ¯ H = RH GH . G c
(5)
C
In the case of a frictionless contact, where the coefficient of friction is either too low or practically zero, the contact modes for a non-conformal contact are reduced to Slip and Separation only [15,13]. The same situation occurs in the case of a sliding contact, where all the contact node pairs slide over each other. Additionally as friction does not affect much the normal pressure distribution [14], the coefficient of friction can be neglected for simplicity, as in the frictionless case. Therefore, considering a sliding wear non-conformal problem, where one body slides over another, the contact modes of the contact node pairs are Slip or Separation only. The boundary conditions implemented for the aforementioned contact modes and the check criteria for any contact mode violations are presented in Appendix B. According to the Multi-Region technique [12], the individual systems of equations are combined for the two (or more) regions and the interface boundary conditions are subsequently applied [13,14]. Considering Eq. (3) for each body H = {A, B}, and after
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
rearranging and applying the pre-scribed boundary conditions of the non-contact boundary, the resulting system of equations can be written as:
tion Wear, which arises from the adhesive forces set up whenever atoms come into intimate contact [17,18]. In delamination wear, the sliding speed is low enough that surface heating can be neglected. The model that describes this mechanism is the well known Archard’s linear wear equation [19], firstly proposed by Holm [20], given as: W H = KH
(6) H where the matrices AH nc and Rnc contain terms of both the H H previous matrices Hnc and Gnc , after the substitution of the prescribed boundary conditions of the non-contact area of each H body H = {A, B}. The column vectors xH nc and ync contain the unknown and the known, respectively, components of tractions and displacements of the non-contact boundary, after the substitution of the pre-scribed boundary values. The submatrix BC1 in Eq. (6) has 2Mc × 8Mc components and contains all the boundary conditions for the displacements and the tangential tractions for the contact node pairs in slip mode and also the constraints for the zero tangential and normal tractions for the pairs in separation mode. The submatrix BC2 has 2Mc × 4Mc components and contains all the boundary conditions for the tractions equilibrium inside the contact-worn area. The column matrices F1 and F2 have 2Mc components each and contain the right-hand sides of the boundary conditions implemented in sub-matrices BC1 and BC2 (Appendix C). The above coupled system, containing the elastostatic problems of both bodies plus the contact problem, is square of size Mt × Mt , where Mt = 2(MA + MB ) + 4Mc . Taking a closer look at Eq. (6) it can be seen that the square matrix is highly sparsed. More than 60% of the total components are zero terms. This becomes significant as wear simulation demands solution of the above system several times, after updating the geometry of the worn body at every step of the solution algorithm. Another advantage of the specific method is that until convergence to be achieved concerning the modes of contact for each node pair, only some components of the submatrix BC1 on the left-hand side of Eq. (6) may change. Moreover as the changes take place only in a small part of the overall system, there is no need to calculate every time the inverse of the total matrix; however, using solution algorithms such as Sherman–Morrison or Woodbury formulas [16] the final solution of the above system can be speeded up several times [12,13].
3. Wear model In the case where the sliding speed remains at low levels and the applied load does not exceed a limit where seizure takes place, the wear that occur is characterized as Delamina-
1121
P S Ho
(7)
where WH is the volumetric material loss of the body in consideration, H = {A, B}, KH the wear coefficient, Ho the hardness of the softer body in contact, P the applied load and SH the sliding distance the body has travelled. For any pair of materials that slide over each other, a specific constant KH can be determined through an experiment. Using a wear model in a computational technique for wear simulation, the parameter “sliding material” is directly introduced in the formulation. However, the model described by Eq. (7) does not provide sufficient details; for that reason instead of the wear volume the Wear Depth can be used for modelling the worn geometry. Moreover, the wear depth is of practical importance in every day engineering as it can be measured directly on a specific mechanical component. In the two-dimensional case, the above model describes the volume of material loss per unit thickness, with respect to the sliding distance, when a specific load P per unit thickness is applied. Therefore, differentiating Eq. (7) with respect to the width of the bearing surface, results to: hH (x, s) = kH p(x, s)sH
(8)
where hH is the wear depth on point x of body H = {A, B} after sH sliding distance, p the pressure distribution developed inside the contact area and kH the modified wear coefficient given by kH = KH /Ho . Both pressure and wear depth are functions of the location and the sliding distance as the bodies lose material during the wear process and the topography of their bearing surfaces changes. Furthermore, in cases as a pin on a rotating disc problem, the travelled distance sH by each point x on the bearing surfaces is also a function of the location of the point. However, this sliding distribution with respect to the location of each point on the bearing surfaces can be evaluated directly by the kinematics of the specific tribosystem. The direction of wear, the wear vector, is taken to coincide with the pressure vector; that is normal to the bearing surface of the body in contact. 4. Solution algorithm The most commonly used technique for simulating wear problems is the Incremental Method [1–7]. It is clear from Eq. (8) that once the pressure distribution inside a contact area is known for a specific sliding distance, the wear depth, with respect to the width of the contact area, can be calculated. However, as the function of the pressure distribution is not known a priori (due to the variation with the sliding distance), the wear depth cannot be calculated directly. To overcome this difficulty the
1122
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
geometry, where the Boundary Element mesh (or the FE mesh) is generated. Then the contact analysis starts with an initial guess of the contact node pairs. The coupled contact problem is solved and checks for contact mode violations are performed. If the initial guess for the contact node pairs that come into contact, for the specific applied normal load, was false, the assumed pairs are updated and the problem is solved again and again until no contact mode violations occur. Then using the resulting pressure distribution inside the contact area, and using Eq. (8), the resulting wear depth is calculated. The geometries of both bodies are updated (in case that both bodies lose material), and the intermediate results are stored. If the accumulated sliding distance reaches the total sliding distance then the simulation is finished. Otherwise a new BE mesh is generated and the above procedure is repeated. Some fine tuning of this method is required in terms of the sliding-time increment that will be used. If the increment is relatively large, then oscillations will occur to the resulting geometries and consequently will not be smooth. Moreover, the validity of the calculated wear depth will not be definite. A trial and error procedure is adopted to calibrate the maximum step that produces smooth geometries and wear depth, and additionally the computational time to be reasonable small. Another way is to use a variable step technique, in order the maximum wear depth for every increment to be less than a specified value [1]; thus, the solution procedure may speed up. 5. Results and discussion 5.1. One body wear formulation
Fig. 2. Solution algorithm for the Incremental Method.
total sliding distance is incremented and the problem is solved for small increments updating every time the geometries of the bodies. The basic algorithm of the Incremental Method using the BEM is illustrated in Fig. 2. The algorithm begins with the initial
In the case where only one of the two bodies in contact loses material, the procedure to simulate a wear experiment is relatively straightforward. To validate the BEM presented here, analytical results using the wear model and the results from the work of P˜odra and Andersson [1] were used. P˜odra and Andersson [1] used the Finite Element Method ANSYS code. In their work 3D formulation was used for simulating the classic problem of a pin on disc wear experiment. Their results were compared, in terms of the wear depth, with the experimental evidences that were used to calibrate the wear model. In the present work, a two-dimensional formulation is employed. For this reason and to allow comparison with the results in [1], the calculated wear depth was transformed to wear volume, using the geometric equation that calculates the volume for a specific depth of a spherical ended pin, given by [21]: WA =
π A 2 A 3R − hA hm m . 3
(9)
For the maximum wear depth hA m calculated using the FEM package, and observed experimentally in the work of P˜odra and Andersson [1], and for a cylindrical ended pin of thickness 1 mm,
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128 Table 1 Parameters used in example number 1 Wear coefficient (Pa−1 ) Total sliding distance (m) Normal load (N/mm) Poisson ratio Radius of the pin (m) Elastic moduli (GPa) Sliding increment (mm)
kA = 1.33 ± 0.54 × 10−13 S = 3.0 P = 50 vA = vB = 0.3 RA = 5 × 10−3 EA = EB = 210 ds = 5, 10, 20, 40, 50
given by [21]: RA − hA m W = (R ) cos RA A 2 − RA − hA 2RA hA m − (hm ) m A
A 2
−1
(10)
for the maximum wear depth hA m calculated using the BEM in the present work. The reason that wear volume is used instead of wear depth to compare the two different formulations (i.e. a 2D and a 3D), is due to the wear law, Eq. (7), as it is not taking into account the geometries of the bodies. Hence, using the same normal applied load, wear coefficient and sliding distance, the wear volume must be exactly the same whatever are the contacting geometries. Each body was modelled using 65 quadratic elements, whereas 69 nodes (34 elements of the 65) were used to model the potential contact-worn area. The parameters of the specific problem are presented in Table 1. The mesh density remains constant throughout the process and consequently the first steps are more crucial since they involve less number of nodes to simulate high
1123
peak contact pressures. As wear progresses, smoother pressures will result and more nodes of the potential contact-worn surfaces will come into contact. Fig. 3 illustrates the initial pressure distribution, before wear, evaluated analytically by Hertzian theory [15] and numerically by the BEM (plane strain). Several trials considering different mesh sizes were carried out to assure the accurate evaluation of the bearing pressure. The mesh size used in the present example (69 nodes for the potential contact-worn area) corresponds to part (c) of Fig. 3, incorporating 9 nodes inside the initial contact area. It is obvious that the developed bearing pressure is evaluated sufficiently accurate, even during the initial step where only 9 nodes are in contact. The same parameters were used in the work of P˜odra and Andersson [1], who calibrate the wear coefficient kA experimentally. In this problem, only the pin is assumed to lose material, as the disc is assumed very hard and its wear can be neglected compared to that of the pin. Moreover, the disc diameter is assumed very large compared to the pin diameter and the position of the pin is far from the center of revolution of the disc. Consequently, as the disc rotates, the difference of the travelled distance between points on the bearing surface of the pin can be neglected and thus all points on the pin, either in contact or not, is assumed to travel exactly the same distance [1]. In case where the travelled distance between points corresponding to the internal and the external radius of the sliding path varies significantly, then the sliding distance distribution must be evaluated from the kinematics of the specific tribosystem, corresponding to every node of the discretized geometry of the body that loses material. The modified wear law, Eq. (8), is applied to all points on the bearing surface of the pin; for points that are in contact and consequently under bearing pressure, the corresponding wear
Fig. 3. Initial contact pressure: (a) #3 nodes in contact; (b) #5 nodes in contact; (c) #9 nodes in contact; (d) #13 nodes in contact (a: width of contact area).
1124
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
Fig. 4. Comparison of the wear volume calculated analytically by the wear model with the wear volume evaluated using the BEM and the Incremental Method, with respect to the sliding distance.
depth is calculated by Eq. (8); however, for points that are not in contact yet, the corresponding wear depth is zero as these points are not under bearing pressure. Therefore, as the pin slides over the disc, the travelled distance by each point on the pin’s bearing surface coincide with the sliding distance of the pin (as a body), while each point on the pin’s bearing surface slides different distance as it comes into contact at different time as the pin loses material. In this paper, the following analysis corresponds to the sliding distance of the pin, as a body. Fig. 4 illustrates the calculated wear volume by the present work, using Eq. (10) in comparison with the analytically calculated wear volume using directly the wear model, Eq. (7) after substituting kH = KH /Ho , and by using the same parameters presented in Table 1. It is apparent that as the sliding increment decreases the wear volume converges. It is obvious that the numerically evaluated wear volume, using the BEM and the incremental form of the wear model, are in good agreement compared with the analytically calculated wear volume. Additionally, in Fig. 5, the maximum wear depth calculated by the proposed formulation is illustrated against the wear depth calculated solving inversely Eq. (10) and considering the wear volume that was evaluated using the wear model, Eq. (7), which
Fig. 6. Comparison of the wear volume evaluated experimentally, by FEM and BEM, with respect to the sliding distance.
is illustrated in Fig. 4. The same behaviour is demonstrated as the sliding increment decreases. Fig. 6 illustrates the wear volume calculated in the present work, compared to the wear volume from [1], by using Eq. (9) for both the numerical and the experimental results that were reported in [1]. The upper and lower limits illustrated in the same figure are the deviation of the wear coefficient k as it was calibrated experimentally by P˜odra and Andersson [1]. In Table 2, the parameters of another numerical example are presented. In this example each body was modelled using 50 quadratic elements in total, whereas 49 nodes (24 elements) were used to compose the potential contact area. Fig. 7 illustrates the final geometry of the worn pin after a total sliding distance of 500 mm, while Fig. 8 illustrates the deviation of the maximum wear depth for various sliding increments. It is obvious from both figures that as the sliding increment decreases the solution converges. However the computational time increases almost linearly as the same calculation processes must be carried out for a greater number of sliding increments. Taking a closer look to Fig. 7, it can be seen that the use of large sliding increments results to the development of an artificial roughness on the bearing surface of the body that wears out. This is due to the incrementation of the sliding distance, as large increments create large craters which produce localized high pressures at their edges after the geometry update at the next step. These localized high pressures create new large craters and the phenomenon continues. As a result geometries with artificial roughness are created. On the other hand, as the sliding increment decreases, these craters are relatively small Table 2 Parameters used in example number 2
Fig. 5. Variation of analytically and numerically calculated maximum wear depth with respect to the sliding distance.
Wear coefficient (Pa−1 ) Total sliding distance (m) Normal load (N/mm) Poisson ratio Radius of the pin (m) Elastic moduli (GPa) Sliding increment (mm)
kA = 1.33 × 10−13 S = 0.5 P = 4.63 vA = vB = 0.3 RA = 10 × 10−3 EA = EB = 70 ds = 2.5, 5, 10, 20
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
1125
Fig. 7. Deviation of the final geometry with respect to the sliding increment: (a) ds = 20 mm; (b) ds = 10 mm; (c) ds = 5 mm; (d) ds = 2.5 mm (a: width of contact-worn area).
and produce a smooth transition from the unworn to the worn geometry. The aforementioned artificially created roughness is also responsible for the deviation of the wear depth and volume, as Figs. 4, 5 and 8 illustrate. 5.2. Two bodies wear formulation In cases where both bodies wear out in comparable time, the problem becomes more complicated as both bearing surfaces change with respect to the sliding distance. In such cases, it may be of interest to simulate the evolution of wear on each body in contact individually. To do so the modified wear law, Eq.
(8), must be used for each body H = {A, B}. Table 3 illustrates the parameters used in the case of two bodies wear formulation. In the specific problem it is assumed that all points on the pin’s bearing surface travel the same distance with all the points on the disc’s bearing surface. This is done for simplicity and to demonstrate the relation between the final wear of each body and its wear coefficient. In practise, in a pin on a rotating disc problem, the sliding distance distribution must be evaluated, corresponding to the changeable location of every node on both bearing surfaces, considering the kinematics of the tribosystem and the contact duration, as it is already explained in Section 5.1. Fig. 9 illustrates the final geometry of the two bodies. Taking a closer look to the specific figure, it can be seen that the maximum wear depth of the pin, body A, is exactly two times the maximum wear depth of the disc, body B. The specific relation between the maximum wear depth of the two bodies results from the input data, Table 3, where the wear coefficient of the pin is exactly two times the coefficient of the disc. In the case where both bodies had the same wear coefTable 3 Parameters for the two-bodies-wear formulation example
Fig. 8. Deviation of the maximum wear depth with respect to the sliding increment.
Wear coefficient (Pa−1 ) Total sliding distance (m) Normal load (N/mm) Poisson ratio Radius of the pin (m) Elastic moduli (GPa) Sliding increment (mm)
kA = 0.887 × 10−13 , kB = 0.443 × 10−13 S = 0.5 P = 4.63 vA = vB = 0.3 RA = 10 × 10−3 EA = EB = 70 ds = 2.5
1126
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
In Fig. 11, the overall maximum wear depth of a wear problem where both bodies lose material, is compared with a problem where only one body lose material, but its wear coefficient is equal to the sum of the two coefficient of the two bodies wear formulation problem. It is obvious that the resulting overall wear depth from the two bodies wear formulation is the same with the one body wear formulation, as the summation of the wear coefficients in the first case is equal to the wear coefficient in the second case. This analysis demonstrates the linear behaviour of the incremented modified wear law, Eq. (8), with respect to the wear coefficients kH . 6. Conclusions Fig. 9. Two bodies wear formulation: Final geometries for k(pin) = 2k(disc) (a: width of contact-worn area).
Fig. 10. Two bodies wear formulation: final geometries for k(pin) = k(disc) (a: width of contact-worn area).
ficients, i.e. kA = kB = 0.665 × 10−13 Pa−1 , the resulting geometries are illustrated in Fig. 10. In this figure the maximum wear depth of the pin is equal to the maximum wear depth of the disc.
In the present work, a two-dimensional application of the Boundary Element Method in wear simulation was presented. The Incremental Method was used where the total sliding distance is divided into small increments. The concluding remarks of the present work are summarized as follows: (1) The stress-strain fields inside the contact area are evaluated fast and accurately using the BEM. Thus, monitoring of the changes in these fields during or after wear, results to a better knowledge of the wear evolution on a specific component. (2) The computational time remains low as only boundary modelling is required in contrast to the FEM. This is a great advantage considering the wear simulation problems, where the geometry is updated after every step. Thus remeshing of the bodies in contact becomes relatively easy and fast. (3) The wear depth is calculated more accurately using the BEM, in comparison with the FEM, as it is directly related to the pressure distribution inside the contact. This is due to the ability of the BEM to evaluate with high accuracy both stresses and displacements unknowns on the boundaries of the bodies in concern. (4) The size of every sliding increment must be optimized in order the final geometries of the bodies to be smooth, in some specified tolerance, and the computational time to remain low. Appendix A The fundamental solutions in Eq. (2) can be obtained by several methods. Here, the fundamental solutions obtained using the Galerkin vector are used [12]. The displacement and traction fundamental solutions for the two-dimensional plain strain condition are given as follows:
1 1 Uij (x , x) = (3 − 4v) ln δij + r,i r,j 8πG(1 − v) r Tij (x , x) =
Fig. 11. Comparison of the maximum wear depth between one and two bodies wear formulation.
∂r [(1 − 2v)δij + 2r,i r,j ] ∂n −(1 − 2v)(r,i nj − r,j ni )
−1 4π(1 − v)r
where Uij (x , x) and Tij (x , x) represents the displacement and the traction, respectively, in the j direction at point x due to a unit
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
point force acting in the i direction at x . G is the shear modulus, nj denotes the components of the outward normal unit vector at the field point x, r is the distance between the field and the force point and δij is the Kronecker delta function. Appendix B The boundary conditions considering separation and slip modes are given in Table 4, where (tcH )t and (tcH )n are components of tangential and normal tractions, respectively, and (uH c )n are components of normal displacements, of each body H = {A, B}. The normal gap between the potential contact node pairs is defined as gapA–B . As shown in Fig. 12 the initial normal gap gapA–B is defined as the normal distance between a node pair, and is calculated as follows [13]: gapA–B = |dy|¯nA nA y + |dx|¯ x
1127
Table 5 Contact mode check criteria Decision . . . assumption
Separation
Contact
Separation Contact
B A−B (uA c )n + (uc )n < gap (tcA or B )n ≥ 0
B A–B (uA c )n + (uc )n ≥ gap (tcA or B )n < 0
Calculating the normal gap in this way, the surfaces of the two bodies will come into contact after the specific gap is closed, either if the nodes are in touch or slide over each other. In order to obtain the correct contact area without any geometrical incompatibility and contact mode violations, the displacements and tractions for every node in the contact zone, must be checked against the criteria presented in Table 5.
(11) Appendix C
¯A for all potential contact node pairs, where n¯ A y and n x are the H components of the average normal unit vector n¯ , of each body H = {A, B}, which is calculated from the outward normal unit vectors of the specific pair of nodes nA , nB by n¯ A =
nA − nB = −n¯ B . |nA − nB |
(12)
Table 4 Boundary constraints Separation
Slip
(tcA )t − (tcB )t = 0 (tcA )n − (tcB )n = 0 (tcA )t = 0 (tcA )n = 0
(tcA )t − (tcB )t = 0 (tcA )n − (tcB )n = 0 (tcA )t = 0 B A–B (uA c )n + (uc )n = gap
Fig. 12. Definition of the average normal unit vector and the initial normal gap.
The sub-system in Eq. (6), corresponding to the sub-matrix BC2 , includes all the boundary constraints for the tractions equilibrium inside the contact area. It can be written as: A ¯t [BC2 ] cB ¯tc ⎤ ⎡ 1 0 0 0 . . . −1 0 0 0 ... ⎥ ⎢ 0 ...⎥ ⎢ 0 1 0 0 . . . 0 −1 0 ⎥ ⎢ ⎢ 0 −1 0 . . . ⎥ = {F2 } ⇒ ⎢ 0 0 1 0 . . . 0 ⎥ ⎥ ⎢ 0 0 −1 . . . ⎥ ⎢0 0 0 1 ... 0 ⎦ ⎣ .. .. .. .. .. .. .. .. . . .. . . . . . . . . . . ⎧ ⎫ ⎧ ⎫ (tcA )n,1 ⎪ ⎪ ⎪ ⎪ ⎪0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (tcA )t,1 ⎪ ⎪ ⎪ ⎪ ⎪ 0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ A ⎪ ⎪ ⎪ ⎪ (t ) ⎪ ⎪ ⎪ ⎪ 0 c n,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ A ⎪ ⎪ ⎪ ⎪ (t ) ⎪ ⎪ ⎪ 0 ⎪ ⎪ c t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ .. ⎪ ⎪ .. ⎪ ⎪ ⎪ ⎨ ⎨ ⎬ .⎬ . × = B ⎪ ⎪ ⎪ (tc )n,1 ⎪ ⎪0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ B ⎪ ⎪ (tc ) ⎪ ⎪0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ t,1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ B ⎪ ⎪ ⎪ ⎪ 0 ⎪ ⎪ ⎪ ⎪ (t ) c ⎪ ⎪ ⎪ n,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ B 0 ⎪ ⎪ ⎪ (tc )t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ . ⎪ ⎩ ⎪ ⎪ .. ⎭ ⎩ .. ⎪ ⎭ . where (tcA )n,1 is the component of the normal traction acting on node 1 of the contact area of body A. This component is related to the normal traction acting on node 1 of the contact area of body B through the boundary condition for the tractions equilibrium (Table 4). Additionally, (tcA )t,1 is the component of the tangential traction acting on node 1 of the contact area of body A and is related to the corresponding tangential traction (tcB )t,1 acting on the coupled node of body B.
1128
G.K. Sfantos, M.H. Aliabadi / Wear 260 (2006) 1119–1128
If the contact node pair 1 is in Slip and the contact node pair 2 is in Separation mode, the sub-system in Eq. (6) corresponding to the sub-matrix BC1 can be written as: ⎧ A⎫ u¯ c ⎪ ⎪ ⎪ ⎪ B⎪ ⎪ A A ⎨ ⎬ [BC1 ]
u¯ c = {F1 } ⇒ [BC1,1 ] A ⎪ ⎪ ¯tc ⎪ ⎪
⎪ ⎩ ¯B ⎪ ⎭
u¯ c u¯ Bc
tc
⎡
1
0
0
0
...
.. .
.. .
.. .
.. .
..
1
0
0
0
.. .
.. .
.. .
.. .
⎢0 0 0 0 ... 0 0 0 0 ⎢ ⎢0 0 0 0 ... 0 0 0 0 ⇒⎢ ⎢ ⎢0 0 0 0 ... 0 0 0 0 ⎣ .
+ [BC1,2 ]
⎧ A ⎫ (uc )n,1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (uA ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ c ⎪ ⎪ A t,1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (u ) ⎪ c n,2 ⎪ ⎪ ⎤⎪ ⎪ ⎪ A ... ⎪ ⎪ ⎪ (uc )t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ...⎥ ⎪ ⎪ . ⎪ ⎥⎪ ⎨ ⎬ . ⎥ . ...⎥ ⎥ (uB ) ⎪ ⎪ c n,1 ⎪ ...⎥⎪ ⎪ ⎪ ⎪ ⎦⎪ ⎪ (uBc )t,1 ⎪ ⎪ ⎪ .. ⎪ ⎪ ⎪ ⎪ . ⎪ ⎪ B ⎪ ⎪ (u ) ⎪ c n,2 ⎪ ⎪ ⎪ ⎪ ⎪ (uB ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ c t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ . ⎩ . ⎪ ⎭ .
⎡
0
0
0
0
... 0
0
0
0
.. .
.. .
.. .
.. .
..
.. .
.. .
.. .
.. .
⎢0 1 0 0 ... 0 0 0 0 ⎢ ⎢0 0 1 0 ... 0 0 0 0 +⎢ ⎢ ⎢0 0 0 1 ... 0 0 0 0 ⎣ .
⎧ A ⎫ (t ) ⎪ ⎪ ⎪ ⎪ cA n,1 ⎪ ⎪ ⎪ ⎪ ⎪ (tc )t,1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ A ⎪ ⎪ (t ) c n,2 ⎪ ⎪ ⎪ ⎤⎪ ⎪ ⎪ A ... ⎪ ⎪ ⎪ (tc )t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ...⎥ ⎪ . ⎪ ⎪ ⎥⎪ ⎨ ⎬ . ⎥ . ...⎥ ⎥ (t B ) ⎪ ⎪ c n,1 ⎪ ...⎥⎪ ⎪ B ⎪ ⎪ ⎦⎪ ⎪ (t ) ⎪ ⎪ ⎪ .. ⎪ ⎪ Bc t,1 ⎪ ⎪ . ⎪ ⎪ ⎪ (tc )n,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ B ⎪ ⎪ (t ) ⎪ ⎪ c t,2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ .. ⎪ ⎭ .
=
⎧ ⎫ gapA−B ⎪ ⎪ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 0 ⎪ ⎪ ⎪ ⎪ ⎨ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩
0 0 .. .
⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
.
¯tc ¯tBc
= {F1 }
References [1] P. P˜odra, S. Andersson, Simulating sliding wear with finite element method, Tribol. Int. 32 (1999) 71–81. [2] P. P˜odra, S. Andersson, Finite element analysis wear simulation of a conical spinning contact considering surface topography, Wear 224 (1999) 13–21. [3] T. MacGinley, J. Monaghan, Modelling the orthogonal machining process using coated cemented carbide cutting tools, J. Mater. Process. Technol. 118 (2001) 293–300. [4] Y.-C. Yen, J. Sohner, B. Lilly, T. Altan, Estimation of tool wear in orthogonal cutting using the finite element analysis, J. Mater. Process. Technol. 146 (2004) 82–91. [5] T.A. Maxian, T.D. Brown, D.R. Pedersen, J.J. Callaghan, A slidingdistance-coupled finite element formulation for polyethylene wear in total hip arthroplasty, J. Biomech. 29 (1996) 687–692. [6] T.A. Maxian, T.D. Brown, D.R. Pedersen, J.J. Callaghan, 3-Dimensional sliding/contact computational simulation of total hip wear, Clin. Orthop. Rel. Res. 333 (1996) 41–50. [7] J.S.-S. Wu, J.-P. Hung, C.-S. Shu, J.-H. Chen, The computer simulation of wear behaviour appearing in total hip prosthesis, Comput. Methods Prog. Biomed. 70 (2003) 81–91. [8] I. Serre, M. Bonnet, R.-M. Pradeilles-Duval, Modelling an abrasive wear experiment by the boundary element method, C.R. Acad. Sci. Paris 329 (2001) 803–808, S´erie II b. [9] Anderson, T., Boundary Elements in two-dimensional Contact and Friction, Diss No 85, Linkoping Institute of Technology, 1982. [10] F. Paris, J.A. Garrido, On the use of discontinuous elements in 2D contact problems, in: Proc. of Boundary Elements VII, Computational Mechanics Publications, Southampton, 1985, 13.27–13.39. [11] G. Karami, R.T. Fenner, A two-dimensional BEM method for thermoelastic body forces contact problems, in: Proc. of Boundary Elements IX, Computational Mechanics Publications, Southampton, 1987, 417–419. [12] M.H. Aliabadi, The Boundary Element Method Applications in Solids and Structures, 2, Wiley, London, 2002. [13] K. Man, M.H. Aliabadi, D.P. Rooke, BEM Frictional Contact Analysis: modelling considerations, Eng. Anal. Bound. Elem. 11 (1993) 77–85. [14] K. Man, M.H. Aliabadi, D.P. Rooke, BEM Frictional Contact Analysis: An Incremental Loading Technique, Comput. Struct. 47 (1993) 93–905. [15] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 2001. [16] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes in Fortran, in: The Art of Scientific Computing, second ed., Cambridge University Press, Cambridge, 1994. [17] S.C. Lim, M.F. Ashby, Wear-mechanism maps, Acta Metall. 35 (1) (1987) 1–24. [18] E. Rabinowicz, Friction and Wear of Materials, second ed., Wiley, New York, 1995. [19] J.F. Archard, Contact and rubbing of flat Surfaces, J. Appl. Phys. 24 (1953) 981–988. [20] R. Holm, Electric Contacts, Gebers, Stockholm, 1946. [21] M.R. Spiegel, J. Liu, Mathematical Handbook of Formulas and Tables, ISE ed., Mc-Graw Hill, New York, 1999.