- Email: [email protected]
Wear simulation with the Winkler surface model
Wear 207 (1997) 79–85
Wear simulation with the Winkler surface model Priit Po˜dra, So¨ren Andersson Machine Elements, Department of Machine Design Royal Institute of Technology, KTH, S-10044 Stockholm, Sweden
Abstract
Recurrent calculations of contact pressure distribution are often the major difficulty in wear simulations of rubbing contacts. The pressure distribution changes with time, as a result of the wear and plastic deformation of the interacting surfaces. To calculate the pressure distribution in a contact, it is often necessary to use advanced numerical methods. One of the most powerful tools is the finite element method (FEM). However, the FEM—as is also the case with most other numerical methods used for contact problems—has the disadvantage of being quite time consuming. Therefore, the advantages and disadvantages of using the simple Winkler surface model in wear simulations have been studied. The Winkler model can be seen as a mattress with many springs with no lateral effects between them. The stiffness of the springs is a critical parameter, which is evaluated for different boundary conditions and cases. It is shown that the sum of the springs’ volumedeformations, called the ‘deformed volume’, can be used to simplify the determination of the contact pressure distribution in a contact. Wear simulations have been made with a Winkler surface model and compared with simulations made previously with an FEM model. The results show a very good agreement. q 1997 Elsevier Science S.A. Keywords:
Wear simulation; Contact pressure; Winkler surface model
1. Introduction
The conventional way to determine the wear properties of tribosystems is to carry out wear experiments. The wear is then determined by measuring the wear depth or volume, either intermittently or on line. The most commonly used laboratory wear test equipment is the pin-on-disc machine. Normally, however, it is difficult to transform the resultsfrom such simple laboratory tests to real operational cases. Thus, field, rig and component tests under realistic running conditions are often preferred, and are also used to a great extent. However, the possibility to understand and predict ‘what causes what’ in real cases is normally quite limited. Therefore, interest in combining field, rig and component tests with simple, well-controlled laboratory tests has increased. Numerical simulations have then been an important tool to improve the ways of comparing test results from different types of wear test, and to improve the use of data from simple laboratory wear tests in simulations of real cases. One advanced and useful numerical tool for surface mechanics problems is the finite element method (FEM). Most commercial FEM software packages, such as ANSYS, provide a large number of options, including material state changes, friction, stress stiffening and large strains. One important disadvantage with FEM, however, is the time and computing-resource consumption. In wear simulations, con-
tact pressure calculations have to be made at each time step. Therefore, it is of great interest to reduce the time necessary for pressure distribution calculations. One possibility is to use the simple Winkler surface model (the model of elastic foundation), with which the contact surfaces are represented by many elastic bars or springs, not connected laterally witheach other.
2. Wear model
The wear model, in which the wear volume rate is proportional to the normal load, was first published by Holm [1] in 1946. The model was based on experimental observation and written in the form QsK
FN H
(1)
where Q (m3 my1) is the wear volume per sliding distance, K is the dimensionless wear coefficient, H (Pa) is the hardness of the softer material and FN (N) is the normal load. The wear coefficient K was initially introduced to provide agreement between theory and experiment. Holm treated it as a constant, representing the number of abraded atoms per atomic encounter.
0043-1648/97/$17.00 q 1997 Elsevier Science S.A. All rights reserved PII S0043-1648(96)07468-6
Journal: WEA (Wear)
Article: 7468
80
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
The wear equation Eq. (1) is often referred to as the ‘Archard wear equation’. Archard refined the theory behind the wear model by assuming the following [2]: c the contact between two surfaces will occur at local interacting asperities; c the real contact area is proportional to the normal load; c a single asperity contact is circular; c for metals, the asperity deformation is plastic; c the system is isothermal. Archard interpreted the non-dimensional wear coefficient K as the probability that an asperity interaction results in the formation of a wear particle. However, that is not the only possible interpretation. It can be treated as the fraction of the total plastically deformed volume [3], or as the fraction of power contributing to the wear formation [4]. Depending on the actual wear mechanism, it can also be taken to reflect the number of load cycles required to form a fragment of an asperity by a fatigue process [5], or to reflect a complex of variables, considering the effects of friction and surfacefilms, for example [6]. Usually, the wear depth h is of more interest than is the wear volume, and Eq. (1) can be rewritten in the form hskps (2) y1 where ksK/H (Pa ) is the dimensional wear coefficient, p (Pa) is the normal contact pressure and s (m) is the sliding distance that the point on the surface slides against the interacting surface. A wear process can be seen as a dynamic process, and the prediction of this process can be seen as an initial value problem. A wear model that corresponds to Eq. (2) may then be described by a differential equation as dh skp (3) ds 3. Winkler surface model
In an elastic contact, the displacement at a point in the contact is affected to different degrees by the pressure and deformation at all points in the contact. This relationship is simplified in the Winkler surface model or elastic foundation model, where the contact surfaces are modelled as a set of elastic bars (Fig. 1) [7]. The shear between them is
neglected, so the contact pressure at a point is only dependent on the actual deformation at that point. The elements are considered as linear elastic bars of length L (m) and with a spring constant or stiffness KN (N my3) in the z direction. We have p KNs uz
(4)
In the general three-dimensional model, the surface elements are regarded as rectangular (2c=2d) bars; in the axisymmetrical model, they are considered as pipes. The contact pressure p over the elements’ contact surface is assumed to be constant. According to the laws of elasticity, the displacement uz of a bar is related to the normal pressure p as EW ps uz L
(5)
where EW (Pa) is a combined elastic modulus of the Winkler contact. Eqs. (4) and (5) give EW KNs L
(6)
3.1. Determination of the elastic foundation stiffness
The stiffness of the bars will be evaluated by comparing the results with known Hertzian solutions of typical nonconformal contacts. We assume that the deformation in the contact is small in comparison with the principal relativeradii of curvature R9 and R0 (m) of the contacting surfaces [8,9]. According to Johnson [7], Eq. (6) can then be reformulated as EU KNsCW bH
(7)
where bH (m) is the calculated shorter semi-axis of the Hertzian contact ellipse and Es[(1yn 21)/E1q(1yn 22)/E2]y1 is the equivalent modulus of elasticity, with n i and Ei (is1, 2) being the Poisson ratios and elastic moduli (Pa) of the contacting materials respectively (CW is a non-dimensional constant). The elastic foundation contact cannot be simultaneously compatible with a Hertzian contact in all aspects. Therefore, the non-dimensional constant CW must be chosen
Fig. 1. Winkler surface model: (a) contacting bodies; (b) schematic view; (c) three-dimensional view.
Journal: WEA (Wear)
Article: 7468
81
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
Fig. 2. CW vs. R9/R0 for different constraints. 3D, general non-conformal contact with the principal relative radii of curvature R9 and R0 of the contacting surfaces; 2D, a plane non-conformal contact.
to be different, by depending on what and how the determined pressure and deformation are going to be used. To determine the usability of the Winkler surface models in wear simulations, the choice of the elastic bar stiffness for different cases will be analyzed. General non-conformal contacts will be considered. In particular, two extreme cases will be studied: (a) the axisymmetrical case, i.e. R9sR0; (b) the plane case, i.e. R0s`. These two cases can be seen as a sphere on a plane, i.e. a point contact, and as a cylinder on a plane, i.e. a line contact. The contact pressure and deformation, determined using the Winkler model, and according to the Hertz theory, are compared. If the contact areas are equal, i.e. (A)Ws(A)H, then the shapes of the contact areas will be different. The Hertzian ellipse is more slender than the Winkler ellipse, except for a case with a sphere against a plane (R9sR0) and a case with a cylinder against a plane (R0s`). For those two cases, the same contact areas are obtained with CWs1.70 (R9sR0) and CWs1.18 (R0s`) respectively. The maximum pressures determined using the Winkler surface model in these two cases are 33% and 18% higher than the corresponding Hertzian pressures. If the length of the shorter semi-axis in a contact is of particular interest, then the lengths determined using the two methods are quite close to each other if CWs1.35. The error will not exceed 7% in either line or point contacts [7]. The maximum pressures are equal, i.e. (pmax)Ws(pmax)H, if CWs0.95 and CWs0.72 for a point contact (R9sR0) and line contact (R0s`) respectively. With the Winkler surface model, the corresponding contact area radius of the point contact will be 16% larger and the semi-width of the line contact will be 18% larger than the values according to the Hertz solution. The compliances are equal, i.e. (d/pmax)Ws(d/pmax)H, in a point contact if CWs0.64. The contact area radius and maximum pressure determined using the Winkler model will
then be 28% larger and 18% lower than the values according to the Hertz solution. To obtain a general view of the different cases, CW is treated as a variable. Fig. 2 shows its dependences on the ratio R9/R0 for different constraints. The differences between the semi-widths a and b, and the maximum contact pressures determined using the Winkler surface model and according to Hertz for different constraints are presented vs. the ratio R9/R0 in Fig. 3. The differences are represented by the ratios aW bW pmax,W Das , Dbs , Dpmaxs aH bH pmax,H
Because the Hertzian contact ellipse is normally more slender than the Winkler ellipse, the long semi-axis, in the general non-conformal contact case, is normally shorter and the short semi-axis is longer than the corresponding Hertzian values (Fig. 3). Fig. 2 indicates that CW can be considered as a constant, independent of the contact geometry in some cases. The errors made are then analyzed for CWs0.72 and 1.35 (Fig. 3(d)). D b increases only marginally but Da decreases and D pmax increases rather sharply with increasing R9/R0. The limits of Db and D pmax for the two-dimensional conditions are also shown in Fig. 3(d). The contact pressure and deformation can be determined using the Hertz equations only for smooth, quite simple contact geometries. For many other contact geometries, the Love–Bjo¨rklund approach [10,11] can be used instead of the Hertz solutions. 3.2. Calculating the pressure distribution
The sum of the loads supported by the elastic elements of the Winkler surface model (Fig. 1) is equal to the applied load FN. We have
Journal: WEA (Wear)
Article: 7468
82
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
Fig. 3. Discrepancies between elastic foundation model and Hertzian results for different constraints: ———, D a; – – –, Db; ———, D pmax; (a) equal contact area; (b) equal pmax; (c) equal (d/pmax); (d) constant CW values. CWs1.35 (bolder lines) and CWs0.72; s, two-dimensional case.
(8)
FNs8Ai pi i
(m2) is the model elements’ surface area. According to Eqs. (4) and (7), the pressure at a point in the contact is given by Ai
psKNuz,
EU bH
where KNsCW
(9)
which, together with Eq. (8), gives FNs8Ai pisKN8Aiuz,i i
i
The element deformations inside the contact area are taken to be positive. The elastic foundation deformed volume VDef (m3) is directly proportional to the applied load, so we have FN VDefs8Aiuz,is K N i
(10)
This relation can be used to determine the pressure distribution by the iterative procedure, as shown in the flow chart in Fig. 4. When the deformed profile geometry is known, the pressures are calculated using Eq. (9). In the three-dimensional models, Ais4cd and, in the axisymmetrical models, Ais4pRid, where Ri is the elements’ average radius and d is the semi-width. In some cases, the iterative procedure can be avoided. The actual surface topography zsf(x, y) can be replaced with a simple geometry zRsf(x, y) (a prism, a cone, a sphere, a cylinder, etc.), which can be the ideal original surface, an envelope of the actual surface or just a geometry that fits the actual surface geometry quite well. Here, it must be ensured
Fig. 4. The iteration approach.
that the element state will not be altered by this replacement. The elements in contact by the actual geometry z(x, y) under a given load must also be in contact by zR(x, y), i.e. no new elements are brought into the contact by the replacement. The space between the actual surface geometry and the replaced simple geometry is called the ‘added volume’ VAdd (m3) and can be both positive and negative (Fig. 5).
Journal: WEA (Wear)
Article: 7468
83
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
4. Results and discussion
This approach was used to simulate the sliding wear of a cylinder on a flat, and that of a sphere on a flat. The friction and the plastic deformations in the contact were omitted. Fig. 5. Possible surface replacement geometries: (a) flat; (b) cone; (c) sphere.
We introduce the term ‘geometrically lost volume’ VGeom (m3) for the volume captured by the replaced surface zRsf(x, y) and the plane of unloaded elastic foundation (Fig. 6). VGeom is equal to the sum of the deformed volume VDef and the ‘worn’ volume. The ‘worn’ volume is the sum of the added volume VAdd and the actual worn volume VWorn. The added volume and the worn volume can be treated in a similar way, because both represent a deviation from an initial surface. We have VGeomsVDefqVAddqVWorn
(11)
VDef is directly proportional to the applied load and is known if the load is known. The sum of VAdd and VWorn is also known, because the difference at each element on the surface between the actual surface and the replaced surface is known. The remaining task is to calculate the penetration of the replaced simple geometry. This penetration dR is a function of VGeom. In the case where the replaced surface coincides with the initial surface and the load has not decreased so that the worn part of the surface will come outside the contact, the determination is quite straightforward. First, the replaced profile penetration is determined as dRsf(VGeom). Next, the element deformations and pressures are calculated. When the pressure distribution is known, a wear increment calculation, based on Eq. (3), is made, i.e. DhiskpiDsi
and the new z value for a surface point i at the simulation step j will be set equal to zi, jszi, jy1qDhi, j
4.1. Cylinder on flat
The radius of the cylinder was Rs3 mm. The contact was modelled using 490 elastic rods, with each rod being 0.6 mm wide. The bodies were both assumed to be made of steel with an elastic modulus of Es210 GPa and Poisson ratio of n s0.3. The normal load was FNs83 350 N my1 and the dimensional wear coefficient was ks2.36=10y18 Pay1. The overall sliding distance was ss22 km and this wasdiscretized into 250 steps, with the sliding increment D ss88 m. The stiffness constant CW was set to 1.35. Fig. 7 shows the wear simulation results with the elastic foundation surface model, and the results of an FEM simulation previously made [12] for the same contact and with the same parameters. 4.2. Sphere on flat
The sphere with a radius of Rs2.78 mm was modelled using 140 axisymmetrical circular-shaped elements, each of which was 2 mm wide. The flat was assumed to be made of steel with an elastic modulus of Es210 GPa and a Poisson ratio of n s0.3; the sphere was assumed to be made of bronze with Es75 GPa and n s0.3. The normal load was FNs20.1 N and the dimensional wear coefficient was ks2.36=10y18 Pay1. The overall sliding distance was again ss22 km and was discretized into 250 solution steps, making the increment Dss88 m. The chosen stiffness constant CW was 1.35. Fig. 8 shows the results from the simulation of the contact modelled as an elastic foundation, and from an FEM simulation of the same contact made previously [12]. The value of the stiffness coefficient CW has a noticeable effect on the initial pressure distribution and the contact area, but has a minor effect on the final worn geometry. The pressure distributions will be quite uniform during the main part of a simulation, so will be fairly independent of the CW value
Fig. 6. Contact of bodies: (a) general view; (b) single element.
Journal: WEA (Wear)
Article: 7468
84
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
Fig. 7. Cylinder-on-flat Winkler model wear simulation: s, comparative FEM results; – – –, Hertz pressure.
Fig. 8. Sphere-on-flat Winkler model wear simulation: s, comparative FEM results; – – –, Hertz pressure.
used. The Winkler model wear simulation results were compared with those previously made by the authors with the commercial FEM program ANSYS 5.0 [12]. The results show very good agreement. However, there are discrepancies during the running-in period. This is expected, because the effect of local deformation is stronger at that stage than later, when the deformation and the pressure distributions are more uniform. The effect of the discretization on the results has also been studied. In the cases simulated, the step lengths of 88 m should not be longer. Longer step lengths will give some numerical instabilities, while shorter lengths do not improve the results. The element sizes have also been varied. Smaller elements than chosen do not improve the results noticeably.
crepancies are during the running-in stage, when the geometrical changes that result from wear are of importance. Later in the process, the pressure distribution will be more uniform. The Winkler surface model can be used to advantage for wear simulations in many applications. The main advantage is the time necessary to determine the deformation and pressure distributions, but the surface can also be divided into smaller elements, which better indicate local phenomena during a wear process.
5. Conclusions
The key parameters of the Winkler or elastic foundation surface model are the size and the stiffness of the elastic elements. The elements used in this work are linear elastic rectangular bars and circular concentric pipes. The choice of element stiffness can achieve pressures and deformations that are comparable with the corresponding Hertzian values—but not simultaneously. The actual case must determine the choice of the stiffness parameter. In wear simulations, the choice is not as critical as expected. The total deformation of the elements in contact is proportional to the applied load for the Winkler surface model. Hence, the penetration can be found by an iterativeprocedure. An analytical method to find the penetration, deformation and pressure distributions is presented, where the actual profile is replaced with the simple profile. The replacement should not change the shape or size of the contact area. Wear simulations with the Winkler model give good agreement with similar FEM simulations. The only noticeable dis-
Appendix A. Nomenclature A a, b c, d CW E FN H h K k KN L p Q R R9, R0 s u VAdd VDef VWorn x, y, z
area (m2) ellipse semi-axes (m) element cross-section semi-lengths (m) stiffness coefficient elastic modulus (Pa) normal load (N; N my1) hardness (Pa) wear depth (m) dimensionless wear coefficient dimensional wear coefficient (Pay1) element normal spring constant (N my3) element length (m) normal pressure (Pa) volume wear per sliding distance (m3 my1) radius (m) principal relative radii of curvature (m) sliding distance (m) deformation (m) added volume (m3) deformed volume (m3) worn volume (m3) Cartesian coordinates (m)
Journal: WEA (Wear)
Article: 7468
85
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
Greek letters
d n
penetration (m) Poisson ratio
Subscripts
H i j
max R W z
1, 2
Hertz solution element encounter simulation step encounter maximum Replaced Winkler model in z direction bodies in contact
References
[1] R. Holm, Electric Uppsala, 1946.
Contacts,
Almqvist & Wiksells Boktryckeri,
[2] J.F. Archard, Wear theory and mechanisms, in M.B. Peterson and W.O. Winer (eds.), Wear Control Handbook, ASME, New York. [3] M.C. Shaw, Dimensional analysis for wear systems, Wear, 43 (1977) 263–266. [4] L. Johansson, Model and numerical algorithm for sliding contact between two elastic half-planes with frictional heat generation and wear, Wear, 160 (1993) 77–93. [5] I.M. Hutchings, Tribology: Friction and Wear of Engineering Materials, 1992. [6] C.N. Rowe, Some aspects of the heat absorption in the function of a boundary lubricant, Trans. ASLE, 9 (1966) 100–110. [7] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1992. [8] S.P. Timoshenko and J.N. Goodier, Theory of Elasticity, 3rd edn., McGraw-Hill, New York, 1970. [9] C. von Horst, Hu¨tte: die Grundlagen der Ingenieurwissenschaften, Springer, Berlin, 1991. [10] A.E.H. Love, Stress produced in a semi-infinite solid by pressure on part of the boundary, Phil. Trans. R. Soc., Ser. A, 228 (1929) 377. [11] S. Bjo¨rklund, A numerical method for real elastic contacts between three dimensional rough surfaces, TRITA-MAE 1995:8, KTH, Stockholm, 1995. [12] P. Po˜dra and S. Andersson, FEA wear simulation of a sliding contact, Proc. OST-94 Symp. on Machine Design, Tallinn, 14–15 April 1994.
Journal: WEA (Wear)
Article: 7468
Wear simulation with the Winkler surface model Priit Po˜dra, So¨ren Andersson Machine Elements, Department of Machine Design Royal Institute of Technology, KTH, S-10044 Stockholm, Sweden
Abstract
Recurrent calculations of contact pressure distribution are often the major difficulty in wear simulations of rubbing contacts. The pressure distribution changes with time, as a result of the wear and plastic deformation of the interacting surfaces. To calculate the pressure distribution in a contact, it is often necessary to use advanced numerical methods. One of the most powerful tools is the finite element method (FEM). However, the FEM—as is also the case with most other numerical methods used for contact problems—has the disadvantage of being quite time consuming. Therefore, the advantages and disadvantages of using the simple Winkler surface model in wear simulations have been studied. The Winkler model can be seen as a mattress with many springs with no lateral effects between them. The stiffness of the springs is a critical parameter, which is evaluated for different boundary conditions and cases. It is shown that the sum of the springs’ volumedeformations, called the ‘deformed volume’, can be used to simplify the determination of the contact pressure distribution in a contact. Wear simulations have been made with a Winkler surface model and compared with simulations made previously with an FEM model. The results show a very good agreement. q 1997 Elsevier Science S.A. Keywords:
Wear simulation; Contact pressure; Winkler surface model
1. Introduction
The conventional way to determine the wear properties of tribosystems is to carry out wear experiments. The wear is then determined by measuring the wear depth or volume, either intermittently or on line. The most commonly used laboratory wear test equipment is the pin-on-disc machine. Normally, however, it is difficult to transform the resultsfrom such simple laboratory tests to real operational cases. Thus, field, rig and component tests under realistic running conditions are often preferred, and are also used to a great extent. However, the possibility to understand and predict ‘what causes what’ in real cases is normally quite limited. Therefore, interest in combining field, rig and component tests with simple, well-controlled laboratory tests has increased. Numerical simulations have then been an important tool to improve the ways of comparing test results from different types of wear test, and to improve the use of data from simple laboratory wear tests in simulations of real cases. One advanced and useful numerical tool for surface mechanics problems is the finite element method (FEM). Most commercial FEM software packages, such as ANSYS, provide a large number of options, including material state changes, friction, stress stiffening and large strains. One important disadvantage with FEM, however, is the time and computing-resource consumption. In wear simulations, con-
tact pressure calculations have to be made at each time step. Therefore, it is of great interest to reduce the time necessary for pressure distribution calculations. One possibility is to use the simple Winkler surface model (the model of elastic foundation), with which the contact surfaces are represented by many elastic bars or springs, not connected laterally witheach other.
2. Wear model
The wear model, in which the wear volume rate is proportional to the normal load, was first published by Holm [1] in 1946. The model was based on experimental observation and written in the form QsK
FN H
(1)
where Q (m3 my1) is the wear volume per sliding distance, K is the dimensionless wear coefficient, H (Pa) is the hardness of the softer material and FN (N) is the normal load. The wear coefficient K was initially introduced to provide agreement between theory and experiment. Holm treated it as a constant, representing the number of abraded atoms per atomic encounter.
0043-1648/97/$17.00 q 1997 Elsevier Science S.A. All rights reserved PII S0043-1648(96)07468-6
Journal: WEA (Wear)
Article: 7468
80
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
The wear equation Eq. (1) is often referred to as the ‘Archard wear equation’. Archard refined the theory behind the wear model by assuming the following [2]: c the contact between two surfaces will occur at local interacting asperities; c the real contact area is proportional to the normal load; c a single asperity contact is circular; c for metals, the asperity deformation is plastic; c the system is isothermal. Archard interpreted the non-dimensional wear coefficient K as the probability that an asperity interaction results in the formation of a wear particle. However, that is not the only possible interpretation. It can be treated as the fraction of the total plastically deformed volume [3], or as the fraction of power contributing to the wear formation [4]. Depending on the actual wear mechanism, it can also be taken to reflect the number of load cycles required to form a fragment of an asperity by a fatigue process [5], or to reflect a complex of variables, considering the effects of friction and surfacefilms, for example [6]. Usually, the wear depth h is of more interest than is the wear volume, and Eq. (1) can be rewritten in the form hskps (2) y1 where ksK/H (Pa ) is the dimensional wear coefficient, p (Pa) is the normal contact pressure and s (m) is the sliding distance that the point on the surface slides against the interacting surface. A wear process can be seen as a dynamic process, and the prediction of this process can be seen as an initial value problem. A wear model that corresponds to Eq. (2) may then be described by a differential equation as dh skp (3) ds 3. Winkler surface model
In an elastic contact, the displacement at a point in the contact is affected to different degrees by the pressure and deformation at all points in the contact. This relationship is simplified in the Winkler surface model or elastic foundation model, where the contact surfaces are modelled as a set of elastic bars (Fig. 1) [7]. The shear between them is
neglected, so the contact pressure at a point is only dependent on the actual deformation at that point. The elements are considered as linear elastic bars of length L (m) and with a spring constant or stiffness KN (N my3) in the z direction. We have p KNs uz
(4)
In the general three-dimensional model, the surface elements are regarded as rectangular (2c=2d) bars; in the axisymmetrical model, they are considered as pipes. The contact pressure p over the elements’ contact surface is assumed to be constant. According to the laws of elasticity, the displacement uz of a bar is related to the normal pressure p as EW ps uz L
(5)
where EW (Pa) is a combined elastic modulus of the Winkler contact. Eqs. (4) and (5) give EW KNs L
(6)
3.1. Determination of the elastic foundation stiffness
The stiffness of the bars will be evaluated by comparing the results with known Hertzian solutions of typical nonconformal contacts. We assume that the deformation in the contact is small in comparison with the principal relativeradii of curvature R9 and R0 (m) of the contacting surfaces [8,9]. According to Johnson [7], Eq. (6) can then be reformulated as EU KNsCW bH
(7)
where bH (m) is the calculated shorter semi-axis of the Hertzian contact ellipse and Es[(1yn 21)/E1q(1yn 22)/E2]y1 is the equivalent modulus of elasticity, with n i and Ei (is1, 2) being the Poisson ratios and elastic moduli (Pa) of the contacting materials respectively (CW is a non-dimensional constant). The elastic foundation contact cannot be simultaneously compatible with a Hertzian contact in all aspects. Therefore, the non-dimensional constant CW must be chosen
Fig. 1. Winkler surface model: (a) contacting bodies; (b) schematic view; (c) three-dimensional view.
Journal: WEA (Wear)
Article: 7468
81
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
Fig. 2. CW vs. R9/R0 for different constraints. 3D, general non-conformal contact with the principal relative radii of curvature R9 and R0 of the contacting surfaces; 2D, a plane non-conformal contact.
to be different, by depending on what and how the determined pressure and deformation are going to be used. To determine the usability of the Winkler surface models in wear simulations, the choice of the elastic bar stiffness for different cases will be analyzed. General non-conformal contacts will be considered. In particular, two extreme cases will be studied: (a) the axisymmetrical case, i.e. R9sR0; (b) the plane case, i.e. R0s`. These two cases can be seen as a sphere on a plane, i.e. a point contact, and as a cylinder on a plane, i.e. a line contact. The contact pressure and deformation, determined using the Winkler model, and according to the Hertz theory, are compared. If the contact areas are equal, i.e. (A)Ws(A)H, then the shapes of the contact areas will be different. The Hertzian ellipse is more slender than the Winkler ellipse, except for a case with a sphere against a plane (R9sR0) and a case with a cylinder against a plane (R0s`). For those two cases, the same contact areas are obtained with CWs1.70 (R9sR0) and CWs1.18 (R0s`) respectively. The maximum pressures determined using the Winkler surface model in these two cases are 33% and 18% higher than the corresponding Hertzian pressures. If the length of the shorter semi-axis in a contact is of particular interest, then the lengths determined using the two methods are quite close to each other if CWs1.35. The error will not exceed 7% in either line or point contacts [7]. The maximum pressures are equal, i.e. (pmax)Ws(pmax)H, if CWs0.95 and CWs0.72 for a point contact (R9sR0) and line contact (R0s`) respectively. With the Winkler surface model, the corresponding contact area radius of the point contact will be 16% larger and the semi-width of the line contact will be 18% larger than the values according to the Hertz solution. The compliances are equal, i.e. (d/pmax)Ws(d/pmax)H, in a point contact if CWs0.64. The contact area radius and maximum pressure determined using the Winkler model will
then be 28% larger and 18% lower than the values according to the Hertz solution. To obtain a general view of the different cases, CW is treated as a variable. Fig. 2 shows its dependences on the ratio R9/R0 for different constraints. The differences between the semi-widths a and b, and the maximum contact pressures determined using the Winkler surface model and according to Hertz for different constraints are presented vs. the ratio R9/R0 in Fig. 3. The differences are represented by the ratios aW bW pmax,W Das , Dbs , Dpmaxs aH bH pmax,H
Because the Hertzian contact ellipse is normally more slender than the Winkler ellipse, the long semi-axis, in the general non-conformal contact case, is normally shorter and the short semi-axis is longer than the corresponding Hertzian values (Fig. 3). Fig. 2 indicates that CW can be considered as a constant, independent of the contact geometry in some cases. The errors made are then analyzed for CWs0.72 and 1.35 (Fig. 3(d)). D b increases only marginally but Da decreases and D pmax increases rather sharply with increasing R9/R0. The limits of Db and D pmax for the two-dimensional conditions are also shown in Fig. 3(d). The contact pressure and deformation can be determined using the Hertz equations only for smooth, quite simple contact geometries. For many other contact geometries, the Love–Bjo¨rklund approach [10,11] can be used instead of the Hertz solutions. 3.2. Calculating the pressure distribution
The sum of the loads supported by the elastic elements of the Winkler surface model (Fig. 1) is equal to the applied load FN. We have
Journal: WEA (Wear)
Article: 7468
82
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
Fig. 3. Discrepancies between elastic foundation model and Hertzian results for different constraints: ———, D a; – – –, Db; ———, D pmax; (a) equal contact area; (b) equal pmax; (c) equal (d/pmax); (d) constant CW values. CWs1.35 (bolder lines) and CWs0.72; s, two-dimensional case.
(8)
FNs8Ai pi i
(m2) is the model elements’ surface area. According to Eqs. (4) and (7), the pressure at a point in the contact is given by Ai
psKNuz,
EU bH
where KNsCW
(9)
which, together with Eq. (8), gives FNs8Ai pisKN8Aiuz,i i
i
The element deformations inside the contact area are taken to be positive. The elastic foundation deformed volume VDef (m3) is directly proportional to the applied load, so we have FN VDefs8Aiuz,is K N i
(10)
This relation can be used to determine the pressure distribution by the iterative procedure, as shown in the flow chart in Fig. 4. When the deformed profile geometry is known, the pressures are calculated using Eq. (9). In the three-dimensional models, Ais4cd and, in the axisymmetrical models, Ais4pRid, where Ri is the elements’ average radius and d is the semi-width. In some cases, the iterative procedure can be avoided. The actual surface topography zsf(x, y) can be replaced with a simple geometry zRsf(x, y) (a prism, a cone, a sphere, a cylinder, etc.), which can be the ideal original surface, an envelope of the actual surface or just a geometry that fits the actual surface geometry quite well. Here, it must be ensured
Fig. 4. The iteration approach.
that the element state will not be altered by this replacement. The elements in contact by the actual geometry z(x, y) under a given load must also be in contact by zR(x, y), i.e. no new elements are brought into the contact by the replacement. The space between the actual surface geometry and the replaced simple geometry is called the ‘added volume’ VAdd (m3) and can be both positive and negative (Fig. 5).
Journal: WEA (Wear)
Article: 7468
83
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
4. Results and discussion
This approach was used to simulate the sliding wear of a cylinder on a flat, and that of a sphere on a flat. The friction and the plastic deformations in the contact were omitted. Fig. 5. Possible surface replacement geometries: (a) flat; (b) cone; (c) sphere.
We introduce the term ‘geometrically lost volume’ VGeom (m3) for the volume captured by the replaced surface zRsf(x, y) and the plane of unloaded elastic foundation (Fig. 6). VGeom is equal to the sum of the deformed volume VDef and the ‘worn’ volume. The ‘worn’ volume is the sum of the added volume VAdd and the actual worn volume VWorn. The added volume and the worn volume can be treated in a similar way, because both represent a deviation from an initial surface. We have VGeomsVDefqVAddqVWorn
(11)
VDef is directly proportional to the applied load and is known if the load is known. The sum of VAdd and VWorn is also known, because the difference at each element on the surface between the actual surface and the replaced surface is known. The remaining task is to calculate the penetration of the replaced simple geometry. This penetration dR is a function of VGeom. In the case where the replaced surface coincides with the initial surface and the load has not decreased so that the worn part of the surface will come outside the contact, the determination is quite straightforward. First, the replaced profile penetration is determined as dRsf(VGeom). Next, the element deformations and pressures are calculated. When the pressure distribution is known, a wear increment calculation, based on Eq. (3), is made, i.e. DhiskpiDsi
and the new z value for a surface point i at the simulation step j will be set equal to zi, jszi, jy1qDhi, j
4.1. Cylinder on flat
The radius of the cylinder was Rs3 mm. The contact was modelled using 490 elastic rods, with each rod being 0.6 mm wide. The bodies were both assumed to be made of steel with an elastic modulus of Es210 GPa and Poisson ratio of n s0.3. The normal load was FNs83 350 N my1 and the dimensional wear coefficient was ks2.36=10y18 Pay1. The overall sliding distance was ss22 km and this wasdiscretized into 250 steps, with the sliding increment D ss88 m. The stiffness constant CW was set to 1.35. Fig. 7 shows the wear simulation results with the elastic foundation surface model, and the results of an FEM simulation previously made [12] for the same contact and with the same parameters. 4.2. Sphere on flat
The sphere with a radius of Rs2.78 mm was modelled using 140 axisymmetrical circular-shaped elements, each of which was 2 mm wide. The flat was assumed to be made of steel with an elastic modulus of Es210 GPa and a Poisson ratio of n s0.3; the sphere was assumed to be made of bronze with Es75 GPa and n s0.3. The normal load was FNs20.1 N and the dimensional wear coefficient was ks2.36=10y18 Pay1. The overall sliding distance was again ss22 km and was discretized into 250 solution steps, making the increment Dss88 m. The chosen stiffness constant CW was 1.35. Fig. 8 shows the results from the simulation of the contact modelled as an elastic foundation, and from an FEM simulation of the same contact made previously [12]. The value of the stiffness coefficient CW has a noticeable effect on the initial pressure distribution and the contact area, but has a minor effect on the final worn geometry. The pressure distributions will be quite uniform during the main part of a simulation, so will be fairly independent of the CW value
Fig. 6. Contact of bodies: (a) general view; (b) single element.
Journal: WEA (Wear)
Article: 7468
84
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
Fig. 7. Cylinder-on-flat Winkler model wear simulation: s, comparative FEM results; – – –, Hertz pressure.
Fig. 8. Sphere-on-flat Winkler model wear simulation: s, comparative FEM results; – – –, Hertz pressure.
used. The Winkler model wear simulation results were compared with those previously made by the authors with the commercial FEM program ANSYS 5.0 [12]. The results show very good agreement. However, there are discrepancies during the running-in period. This is expected, because the effect of local deformation is stronger at that stage than later, when the deformation and the pressure distributions are more uniform. The effect of the discretization on the results has also been studied. In the cases simulated, the step lengths of 88 m should not be longer. Longer step lengths will give some numerical instabilities, while shorter lengths do not improve the results. The element sizes have also been varied. Smaller elements than chosen do not improve the results noticeably.
crepancies are during the running-in stage, when the geometrical changes that result from wear are of importance. Later in the process, the pressure distribution will be more uniform. The Winkler surface model can be used to advantage for wear simulations in many applications. The main advantage is the time necessary to determine the deformation and pressure distributions, but the surface can also be divided into smaller elements, which better indicate local phenomena during a wear process.
5. Conclusions
The key parameters of the Winkler or elastic foundation surface model are the size and the stiffness of the elastic elements. The elements used in this work are linear elastic rectangular bars and circular concentric pipes. The choice of element stiffness can achieve pressures and deformations that are comparable with the corresponding Hertzian values—but not simultaneously. The actual case must determine the choice of the stiffness parameter. In wear simulations, the choice is not as critical as expected. The total deformation of the elements in contact is proportional to the applied load for the Winkler surface model. Hence, the penetration can be found by an iterativeprocedure. An analytical method to find the penetration, deformation and pressure distributions is presented, where the actual profile is replaced with the simple profile. The replacement should not change the shape or size of the contact area. Wear simulations with the Winkler model give good agreement with similar FEM simulations. The only noticeable dis-
Appendix A. Nomenclature A a, b c, d CW E FN H h K k KN L p Q R R9, R0 s u VAdd VDef VWorn x, y, z
area (m2) ellipse semi-axes (m) element cross-section semi-lengths (m) stiffness coefficient elastic modulus (Pa) normal load (N; N my1) hardness (Pa) wear depth (m) dimensionless wear coefficient dimensional wear coefficient (Pay1) element normal spring constant (N my3) element length (m) normal pressure (Pa) volume wear per sliding distance (m3 my1) radius (m) principal relative radii of curvature (m) sliding distance (m) deformation (m) added volume (m3) deformed volume (m3) worn volume (m3) Cartesian coordinates (m)
Journal: WEA (Wear)
Article: 7468
85
P. Po˜dra, S. Andersson / Wear 207 (1997) 79–85
Greek letters
d n
penetration (m) Poisson ratio
Subscripts
H i j
max R W z
1, 2
Hertz solution element encounter simulation step encounter maximum Replaced Winkler model in z direction bodies in contact
References
[1] R. Holm, Electric Uppsala, 1946.
Contacts,
Almqvist & Wiksells Boktryckeri,
[2] J.F. Archard, Wear theory and mechanisms, in M.B. Peterson and W.O. Winer (eds.), Wear Control Handbook, ASME, New York. [3] M.C. Shaw, Dimensional analysis for wear systems, Wear, 43 (1977) 263–266. [4] L. Johansson, Model and numerical algorithm for sliding contact between two elastic half-planes with frictional heat generation and wear, Wear, 160 (1993) 77–93. [5] I.M. Hutchings, Tribology: Friction and Wear of Engineering Materials, 1992. [6] C.N. Rowe, Some aspects of the heat absorption in the function of a boundary lubricant, Trans. ASLE, 9 (1966) 100–110. [7] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1992. [8] S.P. Timoshenko and J.N. Goodier, Theory of Elasticity, 3rd edn., McGraw-Hill, New York, 1970. [9] C. von Horst, Hu¨tte: die Grundlagen der Ingenieurwissenschaften, Springer, Berlin, 1991. [10] A.E.H. Love, Stress produced in a semi-infinite solid by pressure on part of the boundary, Phil. Trans. R. Soc., Ser. A, 228 (1929) 377. [11] S. Bjo¨rklund, A numerical method for real elastic contacts between three dimensional rough surfaces, TRITA-MAE 1995:8, KTH, Stockholm, 1995. [12] P. Po˜dra and S. Andersson, FEA wear simulation of a sliding contact, Proc. OST-94 Symp. on Machine Design, Tallinn, 14–15 April 1994.
Journal: WEA (Wear)
Article: 7468