of linearly reciprocating block-on-flat sliding test
André Rudnytskyj
Mechanical Engineering, master's level (120 credits) 2018
Luleå University of Technology
Department of Engineering Sciences and Mathematics
Andr´e Rudnytskyj
Master programme in Tribology of Surfaces and Interfaces - TRIBOS 4th ed.
Enrollment Semester Autumn 2016.
Department of Engineering Sciences and Mathematics Lule˚a University of Technology
©Lule˚a University of Technology, 2018. This document is freely available at www.ltu.se
ment of Engineering Sciences and Mathematics of Lule˚a University of Technology (LTU), in Sweden. It was made possible through the Joint Erasmus Mundus Mas- ter Course (EMMC) TRIBOS (Tribology of Surface and Interfaces), of which I took part in its 4th generation between the years of 2016 to 2018.
I would like to thank the Education, Audiovisual and Culture Executive Agency (EACEA) of the European Union and the TRIBOS programme coordinators, pro- fessors, tutors, lecturers, and everyone involved in the programme from the Uni- versity of Leeds (UK), the University of Ljubljana (Slovenia), the University of Coimbra (Portugal), and Lule˚a University of Technology (LTU) for their support inside and outside the classroom. I also thank the people directly involved in the development of the thesis for their helpful support and discussions throughout the time of this work.
Andr´e Rudnytskyj Lule˚a, June 2018
with engineering problems, ultimately saving time and resources. In this work, a model problem and methodology is developed to deal with a common situation found in experiments in tribology, namely a linearly reciprocating block-on-flat dry sliding contact. The modelling and simulation of such case would allow a better understanding of the contact pressure distribution, wear and geometry evolution of the block as it wears out during a test. Initially, the introduction and motivation for this work is presented, followed by a presentation of relevant scientific topics related to this work. Wear modelling of published studies are reviewed next, along with studies available in the literature and the goals for this thesis.
The fourth section refers to the methodology used and the built-up of the model problem. In this work the Finite Element Method and Archard’s wear model through COMSOL Multiphysics® and MATLAB® are used to study the proposed contact problem. The construction of the model problem is detailed and the procedure for wear, geometry update and long term predictions, is presented inspired by the literature reviewed.
Finally, the results are presented and discussed; wear increment and new ge- ometries evolution are presented in the figures, followed by pressure profile evo- lution at selected times. The final geometry is also compared for different time steps. At last, conclusions and recommendations for future work are stated.
1.1 Motivation . . . 1
2 Theory & Literature Review 2 2.1 Tribology . . . 2
2.2 Friction . . . 3
2.3 Wear . . . 4
2.4 Factors affecting friction and wear . . . 5
2.4.1 Contact area . . . 5
2.4.2 Load . . . 6
2.4.3 Sliding velocity & Temperature . . . 6
2.4.4 Running-in & Transfer-film . . . 7
2.5 Wear modelling . . . 9
2.5.1 Contact modelling . . . 9
2.5.2 Finite Element Method . . . 11
2.5.3 Archard’s wear equation . . . 12
2.5.4 Wear modelling in the literature . . . 13
3 Research Gap 18 4 Model & Methods 19 4.1 Software . . . 19
4.2 Model problem . . . 19
4.2.1 Wear & Geometry update . . . 24
4.3 Long term predictions . . . 26
4.4 Methodology overview . . . 27
5 Results & Discussion 29 5.1 Contact pressure . . . 29
5.2 Wear & Geometry . . . 31
5.2.1 Initial simulations. . . 32
5.2.2 Finer time steps. . . 34
5.2.3 Full-scale simulation . . . 39
6 Conclusions & Future works 43
cation regimes along with scheme of physical cases; I: Boundary, II:
mixed, III: hydrodynamic, IV: EHL (non-conformal contact only).
Adapted from [33]. . . 2 2 Representation of the real contact area of rough surfaces in contact,
composed by a number of micro-contacts (in red) defined by ran- dom distribution of numerous small contact points, which prevents interlocking or meshing. Side view in square and contact surface view in circle. Adapted from [86]. . . 5 3 Schematics illustrating evolution of contact between a solid lubri-
cant and a rough hard counterpart: from left to right - before sliding, gradual filling of the solid lubricant in the grooves of the counter-surface, buildup of the solid lubricant leading to self lu- brication. Note that there exists a change in roughness of both surfaces. Adapted from [81, 51]. . . 8 4 Point mass supported by spring system (a); illustration of the La-
grange Multiplier Method (b); penalty spring due to penalty term (c). Adapted from [93]. . . 9 5 A plate under uniaxial tension with a hole; finite element (FE) mesh
of the model (middle) with finer (smaller) elements in regions where the fields and errors are expected to be high; typical visualization of a simulation. Adapted from [25]. . . 11 6 Wear occurs when a fragment of material detaches from an asperity
during contact. The shape is approximated as a hemisphere with radius a. Adapted from [86]. . . . 12 7 Geometric model of disc brake assembly on which wear of the pad
was modelled and numerically studied; Surface profile of brake pad for wear simulation. From [84]. . . 15 8 Left: illustration of block on plate experiment (from [28]); right:
modelling of experiment - refer to text for numbered explanation. . 21 9 FE mesh for the counter surface showing finer refinement towards
region where contact with the block occurs. . . 21 10 Finalized initial FE mesh showing zoomed-in details of finer refine-
ment on destination boundary (block). . . 23 11 High peaks and gradients of contact pressure in the block’s edge
region, leading to the set-up of finer FE mesh. . . 24 12 Detail on the construction of the block’s contact boundary; higher
density of points in the region of high pressure peaks and gradients. 25 13 Flowchart of wear computation methodology. . . 28 14 Pressure profile in elastic contact for different indenters (sphere, flat
die and cone). From [96]. . . 29 15 Reference entities for non-dimensionalization of results. . . 29 16 Contact pressure profile (top) normalized by the maximum pres-
sure of the frictionless case, pr,f ric (left) and x-component of the displacement field (bottom); frictional case with µ = 0.3 - chosen for visualization purposes only (right). . . 30
cient of friction µ changes. . . . 31
19 Wear, ∆tg = 2s (“cycle-update”), A = 1000, 1st cycle. . . . 32
20 Wear, ∆tg = 2s (“cycle-update”), A = 1000, 2nd, 3rd and 4th cycles. 33 21 Wear, ∆tg = 2s (“cycle-update”), A = 500, 4th cycle. . . . 34
22 Wear, ∆tg = 0.05s, A = 1000, at different moments. . . . 35
23 Maximum dynamic pressure during simulations of 4 cycles; “*” refers to inversion in the direction of motion when associated with a step “s”, or start of a new cycle when associated with a cycle “c”. 36 24 Dynamic pressure, ∆tg = 0.05s (named dt on the title), A = 1000, in the beginning and in the end of the same stroke moving to the right direction. . . 37
25 Block’s edge final geometry after 4 cycles with A = 1000. . . . 38
26 Block’s edge final geometry after 4 cycles with A = 500. . . . 38
27 Block’s edge final geometry after 2 cycles with A = 1000. . . . 39
28 Maximum dynamic pressure during simulations of first 4 cycles with A = 4000; “*” refers to inversion in the direction of motion when associated with a step “s”, or start of a new cycle when associated with a cycle “c”. . . 40
29 Maximum dynamic pressure during simulations of 4th to 9th cycles with A = 4000. . . . 40
30 Step showing non-smooth wear increment during the 6th cycle with A= 4000. . . 41
31 Wear, ∆tg = 0.05s, A = 4000, at different moments. . . . 42
1 Introduction
Essentially all contacts of solid matter involve friction and wear, from modular junctions of hip replacement, to bearings in turbines of hydropower plants, and NASA’s Mars rover Curiosity’s wheels. Wear is a crucial factor that determines the service life-span and reliability of engineering systems, having enormous eco- nomical impact. In fact, a few decades ago, mechanical wear was estimated to be responsible for a loss of as high as 6 percent of the annual U.S. gross domestic product [22, 69].
In such scenario, scientific and empirical efforts on expanding the understand- ing of wear processes have been made throughout the years, by mathematically developing relationships with relevant parameters in the tribological system. Such challenging and bold endeavours have been exposed by Meng and Ludema [56], whose work reviewed hundreds of wear and friction models from the literature. Al- though for many specific situations, models and theory have been detailed fairly well, there is evidently no universal model applicable to all cases.
Among varied approaches to model events of the natural world, numerical mod- elling comes up as a useful, fast, cheap, improvable complement to experiments, becoming interesting tools to provide scientists with better understanding of a system’s behaviour and events, which can result in accurate predictions providing the model is validated. A general trend towards increased use and development of numerical modelling in several instances of engineering, allow its use also in tribology.
1.1 Motivation
While modelling and simulation of single-part bodies and mechanical evaluation are simulated fairly accurate, assemblies containing physical interfaces and mate- rial loss can raise the difficulty of the task. Due to the complex, multi-disciplinary, and stochastic nature of tribological processes, i.e., surface topography, frictional heat, debris, chemical reactions, atomic interactions and so on, wear modelling is a tricky phenomenon to model, to say the least.
Despite the complexity of computational tribology yielding research in the de- velopment of numerical methods and robust algorithms themselves for implemen- tation of simple contact and friction, established computational tools are readily available to model frictional contact problems. Such tools can be sufficient for study purposes of practical cases, bearing in mind the limitations of the modelling the tool can provide.
Oscillatory contacts are commonly encountered in engineering applications and also in tribology experiments; coefficient of friction and wear coefficient are often estimated in experiments with tribometers such as the pin-on-disc or block-on-flat (also named flat-ended pin-on-disc, [39]) tests, where a sample is pressed against a counter-face and slide with a reciprocating (repetitive up-and-down or back- and-forth) motion. A numerical model that would capture the phenomena and evolution of events taking place in such experiments would serve as a doorway to a better comprehension of the experiment and tribology in general, besides improving design capabilities, overcoming limitations of experiments and saving resources.
2 Theory & Literature Review
It is essential to have a good understanding of the phenomena that may occur in the tribological contact of a linearly reciprocating sliding system, not only for the comprehension of the events in experiments and real applications, but also for the ability to critically analyze the validity of simulations.
In this section, tribology and its relevant topics for this work are presented, followed by a brief introduction to computational tribology and a review on wear modelling by numerical means.
2.1 Tribology
Tribology is the “science and technology of interacting surfaces in relative motion and of the practices related thereto” or, in other words, the study on lubrication, friction, and wear of moving or stationary parts in physical contact (essentially rolling, sliding, normal approach or separation of surfaces) [33, 86]. Tribology comprises not only mechanical phenomena, but several topics in materials science and chemistry, making it a greatly interdisciplinary field.
The lubricant, with the general purpose of - but not limited to - reducing energy loss due to friction and reducing wear is a key part of the tribosystem. Depending on the application, the lubricant may be solid or fluid, and the lubrication regime is classified as: boundary, mixed (or partial), elastohydrodynamic (EHL), and full film lubrication (or hydrodynamic). More details of the regimes can be found, for example, in [86].
The German engineer Richard Stribeck, presented a useful concept to visualize the behaviour of friction in each lubrication regime, in what became known as the Stribeck curve [33]. The curve (Figure 1) present the coefficient of friction under steady state conditions as a function of the Hersey number: Hs =ηω/p, where η is the absolute viscosity, ω is rotational speed, and p is pressure. The higher is Hs, the thicker is the film thickness between the surfaces.
Coefficientoffriction
Hersey number or film parameter Absence of
boundary lubricants
II III I
IV
Figure 1: Stribeck curve showing the behaviour of friction for different lubrication regimes along with scheme of physical cases; I: Boundary, II: mixed, III: hydrodynamic, IV: EHL (non-conformal contact only). Adapted from [33].
In the case of rough surfaces, the Stribeck curve is commonly characterized by λinstead of Hs[86,54], where λ is a specific film thickness with respect to surface roughness of the bodies; the film thickness can be determined from calculations of EHL film in terms of dimensionless entities of speed, materials, load and geometry, famously studied by Hamrock & Downson [33].
In this work, the focus is on boundary lubrication (or boundary film lubrica- tion), which is the scenario of many practical applications such as metal-working tools or polymer bearings in hydropower applications. This regime has as main features [36]:
• Surfaces of the bodies are in contact at microasperities, hence the load is mainly transmitted by mechanical contact resulting in stresses and deforma- tion of the body(ies);
• Tribology is defined by physical properties of the bodies, and possibly by chemical properties of thin surface films of molecular proportions (1-100nm).
2.2 Friction
Friction can be defined as the resistance encountered when one body moves tan- gentially over another with which it is in contact, being the major cause of wear [86] and energy dissipation. It is governed by processes occurring mainly in the thin surface layers of the contacting bodies and thus not an intrinsic material property, but rather of the tribological system. Therefore, when studying and measuring friction, it is important to keep in mind its stochastic nature.
Boundary lubrication regime can be described by means of the coefficient of friction, µ, defined as
µ= Ft
F (1)
where Ft is the tangential frictional force and F is the normal load to the sur- face. The nature of Ft essentially results from two major processes1 in common materials:
• Adhesion, or, to be more precise, the shearing of adhesive junctions, which are created at the interface and constitute the real area of contact;
• Deformation of the softer material caused by penetration of hard asperities of the mating surface, i.e., ploughing of a groove, which can be partly elastic and partly plastic, depending on the materials and roughnesses of the bodies involved [48].
The idea of adhesion and deformation being non-interacting components of a two-term model of friction is found in the literature as a convention; other authors include the separation and formation of transfer film and wear debris as a third factor [60, 59].
The contribution of each process to the total frictional force is a function of several factors; the shearing of adhesive junctions depends, for instance, on the
1Other processes are separately considered and may be more or less important in specific cases, such as scratching, slip, or hysteresis losses [30].
nature of the contact, the chemistry, the stress state and so on, whereas the deformation process is closely related to the mechanical properties of the surfaces, the sliding and environmental conditions. Adhesion can be described in terms of surface energy, the junction growth occurs in a limited way, but, both in polymers and metals, rupture tends to occur within the bulk of the polymer rather than at the interface.
The general understanding of friction between unlubricated solids consists of three basic elements detailed by Tabor [88]: true area of contact; type and strength of bond formed at the interface; and the manner in which the material in and around the contacting regions is sheared and ruptured. Hence, such model assumes that frictional effects occur in a microscale.
2.3 Wear
Wear can be generally defined as the progressive damage to contact areas in relative motion, causing gradual loss of material. Evidently, the factors described in the previous subsections that affect friction will also have an effect on wear. There is a great diversity of mechanisms and conditions related to wear which make a rigorous classification inviable.
In abrasive wear, a harder material/body removes material from a softer body, which can occur by cutting, microfracture on brittle abraded materials, pull-out of individual grains, or repeated deformations by subsequent grits [86]. Cutting refers to the classical model, resulting in wear grooves, micro-machining, or scratching.
The abrasive particles can be present in the microstructure, enter as contaminants, or generated as wear debris, the abrasive action can be in terms of ploughing, i.e., plastic grooving without material removal, or in terms of cutting. The process can be thought of consisting of three stages [48]: mechanical deformation and stresses at a contact area; action of the frictional forces; rupture of material.
Adhesive wear occurs during sliding due to the continuous sticking and tearing of the adhesive junctions, thus it is related to interfacial wear and friction due to shearing of the junctions. Abrasive wear, on the other hand, is related to the ploughing phenomenon. This overview follows the concept of the two-term friction model, i.e., that the energy dissipation of the two frictional terms results in two general modes of damage and consequent wear: interfacial, and cohesive [12]. The balance between interfacial bonding and cohesive forces of the material determine how the fracture will occur.
Repeated stressing during reciprocating can lead to subsurface crack initiation which may propagate parallel to the surface due to deformation. This event results in surface fatigue wear, a common type of wear in polymer tribological systems.
Although it is regarded as a mechanical process of failure, it can also be treated as
“a thermally activated failure process, where failure occurs as a result of thermal breakdown of the chemical bonds and the magnitude of the applied stress reduces the activation energy barrier for bond rupture ” [48]. Evidently, the intensity of failure depends on the materials involved; Smurugov et al. [83] studied how the wear rate for PTFE varies with temperature at different loads, justifying the results in terms of changes in the energy barrier of thermodegradation of the polymer.
Several other types and classification of wear exist, which may involve one or more mechanisms and present resemblances with the previous discussed types.
Stachowiak and Batchelor [86] discuss several aspects of erosive wear, where a surface is damaged by the impact of solid or liquid particles of solid against the surface of an object; cavitation wear, corrosive wear, oxidative wear, fretting wear and other minor mechanisms also exist [39].
2.4 Factors affecting friction and wear
The tribological system is sensitive to numerous factors, which range from material properties to operation conditions. The contact at the interface, its bonds and rupture are all affected by the conditions, consequently affecting friction as a whole.
No single wear mechanism can completely describe wear in a tribological system;
rather, different mechanisms may take place at different times and may become more or less important, slowly or abruptly, as the sliding conditions change. The main important factors affecting friction and wear are considered in the following subsections.
2.4.1 Contact area
It is well known that the real contact area between solid bodies is generally only a small fraction (Figure 2) of the nominal contact area [21] which determine the nature of friction and wear.
Figure 2: Representation of the real contact area of rough surfaces in contact, composed by a number of micro-contacts (in red) defined by random distribution of numerous small contact points, which prevents interlocking or meshing. Side view in square and contact surface view in circle. Adapted from [86].
One of the classical friction laws2 refers to the independence of friction to the size of the nominal contact area3. Indeed, with a large nominal contact area, the micro-contacts (or asperities) lie further apart, hence the sum of these micro- contacts, which is the real contact area, Ar, is, in this sense, independent of the nominal contact area. Greenwood and Williamson famously developed a model
2It is known that ancient civilizations had a substantial understanding of friction and tribology through practical problems associated thereto; empirical laws of sliding friction from what might be termed scientific studies of tribology, often known as Amontons’ Laws, are a re-discovery of work done by Leonardo Da Vinci.
3The other 2 classical friction laws are: proportionality of friction to normal load, and inde- pendence of sliding velocity, the latter added by Coulomb. They have no physical basis and are empirical observations [65].
for multi-asperity contact area by joining Hertzian theory of elastic contact and probability theory, concluding that under certain conditions the ratio of contact area and load can be found to be practically constant [31].
According to Hutchings and Shipway [39], polymers under sufficiently high loads will result in Ar being equal to the nominal contact area, allowing a direct comparison between a measured coefficient of friction and measured shear yield strength of bulk samples of the material. There is no general rule as to the behaviour of friction as a function of roughness. In a polymer-on-metal contact, the metal counter-surface roughness, shape and orientation of the asperities may play an important role in both friction and wear. If a metal counter-surface is very smooth, adhesion forces may become strong, whereas higher surface roughness may result in more abrasion [1].
In this sense, one may realize that the coefficient measured in experiments is only a statistical average, where it is assumed a uniform distribution of shear resistance over all the contact spots, which is not the case for real interfaces.
2.4.2 Load
Another of the classical friction laws refers to the proportionality of friction to the load. For materials with low Young’s modulus to hardness values (E/H), such as common polymers, micro-contacts will largely deform elastically, unless very rough surfaces are involved ([92, 39]). In such case, Ar increases by less than the proportional amount with the normal load, because of increasing com- pression stress in the micro-contacts, but not in the number of asperity contacts;
from Hertzian contact theory [43, 31], any elastic contact in which the number of contacts remains constant under a compressive load F results in
Ar ∝ F2/3 (2)
and, in the case of perfectly plastic behaviour of asperities,
Ar ∝ F. (3)
In the case of a thin polymer film between rigid surfaces, the coefficient of friction decreases with increasing normal load, which is observed in practice for polymers which form transfer films [75]; in fact, it is expected that µ ∝ F−1/3, proportion found in fair agreement in practice for high loads or smooth surfaces [39]. In metals, on the other hand, micro-contacts’ deformation is mostly plas- tic which implies in an increase in the number of micro-contacts as the load is increased (the yield pressure for each asperity contact becoming a material con- stant). Since adhesion is also proportional to Ar, the total friction ends up being roughly proportional to the load [88]. Further discussion can be found in [60]. In terms of wear, the increasing load will result in higher contact pressures and shear stresses, which eventually leads to more mechanical damage.
2.4.3 Sliding velocity & Temperature
With regards to sliding velocity, ploughing and adhesion causing friction in metals are generally independent of velocity; however, when high velocities are involved, melting of asperities may occur which lead to a dependence of the coefficient of
friction on velocity, decreasing with increasing sliding speeds. Furthermore, slid- ing contact in general is “distributed over a lesser number of larger contact areas rather than a large number of contact points” [86], which may be explained by greater separation between opposing surfaces during sliding, depending on the case. Polymers may display a strong visco-elastic behaviour, resulting in an in- crease of E at higher deformation rate [76]; hence, summits penetrate less in the polymer decreasing ploughing friction and real contact area. It has been shown, however, that many polymers exhibit higher wear at higher sliding speeds [75].
Temperature rise normally occurs as the sliding speed increases, since the dis- sipation of the frictional energy in the rubbing is mostly in the form of heat. In metals, limited temperature changes (< 150°C) do not significantly change their material properties, but higher temperatures may cause changes such as decrease in hardness and melting of asperities, increase in rate of oxidation, or even phase transformation, depending on the dissipation. In polymers, mechanical properties may decline fast with temperature rise, as it transits from a “glassy state” (high strength stiffness) to a “rubbery state” (low strength and stiffness); such softening leads Ar to approach the nominal contact area. Hence, on one hand, tempera- ture increase might increase the contact area, but on the other it decreases shear strength. In some thermoplastics, over a critical sliding speed, the temperature at the interface may cause surface melting and thermal softening, leading to slight decrease in wear rate ([1]), while in other materials, the wear rate may increase due to higher adhesion. Evidently, the length scale being study dictates the rel- evance of the events to the general scenario; on the fine scale, for instance, there will likely be a transformation of a contacting layer, whereas in a bigger scale, thermal stresses may relevant in the general load configuration of the system.
Temperature in tribological contacts is a complicated topic, to say the least.
Peaks of temperature called flash temperature are considered to occur in the in- stantaneous true contact area due to adhesive friction in a very short period of time. The many available models to calculate the contact temperature consider different geometrical, physical, and dynamic assumptions, and can vary several hundreds of °C among each other under the same input parameters, as studied by [44]. Measurements encounter the difficulty of not only the contact location being enclosed between two surfaces, but also the temperature varying spot-to-spot in the asperity micro-contacts within the same contact.
2.4.4 Running-in & Transfer-film
The initial stages of a contact between two bodies are commonly characterized by simultaneous transitional processes occurring on different scales within the in- terface; such processes depend on how the energy in the tribosystem is dissipated and are commonly related to friction, wear, temperature and, in some cases, third- body distribution, crystallographic orientation and state of work-hardening. The net result of all these processes is named running-in, or breaking-in, but the term wearing-in refers to when the wear rate during the early life of a machine compo- nent is changed due to geometrical changes. Hence, running-in, besides possibly involving wearing-in, also includes changes in friction that do not necessarily take place over the same period of time as wearing-in; further discussion on this issue can be found in [10].
The sliding process can generate third-body material primarily from the softer bodies, which can be pushed out of the contact interface (wear debris), or be transferred to the harder counterpart and adhered to the surface. In a contact involving polymers (but not exclusively), a phenomenon of material transfer, usu- ally referred to as transfer film formation (or third-body layer), may occur at the interface of the contact affecting the tribological behaviour of the pair. As polymer particles are worn out, depending on the static electric attraction, me- chanical embodiment and sliding conditions of some materials, they form strong adherence to the counter-surface at favourable points due to van der Waals forces, Coulomb electrostatic forces, and chemical bonding [26]. This transferred film increases with the steady rubbing motion resulting in changes in friction and wear by means of shielding of the soft polymer surface from the hard metal asperities, which also play a role in the process by interlocking polymer fragments [41,6,49].
Hence, the film formation is governed by the material composition, counter-surface characteristics and sliding conditions [95].
The capacity to transfer microscopic amounts of material to the mating sur- face creating a film that provides lubrication and reduces friction is known as self-lubrication (Figure 3). The stages of an ideal tribofilm formation is presented by [26,2]: initially there is slight chemical modification in a non-equilibrium zone;
then, an elastoplastic tribofilm is formed, which may involve wear; finally, a sus- tainable tribolfilm is continuously removed and formed during sliding achieving self-lubrication. The behaviour of the transfer is reviewed by [12], dealing with whether the transfer is isothermal or adiabatic, and whether there might be soften- ing, chemical change, morphology change (smooth or lumpy layer), or degradation in the transfer layer.
Figure 3: Schematics illustrating evolution of contact between a solid lubricant and a rough hard counterpart: from left to right - before sliding, gradual filling of the solid lubricant in the grooves of the counter-surface, buildup of the solid lubricant leading to self lubrication. Note that there exists a change in roughness of both surfaces. Adapted from [81,51].
In a sense, the phenomenon of transfer film formation can be seen as a wear process, desirable in some applications [7]. The interfacial processes due to ad- hesion and related to the transfer film are temperature rise, polymer softening, melting, deformation, diffusion, and degradation. It is important to highlight that these many phenomena may occur simultaneously at the contact.
Remark 1. Other factors such as size of wear particles, humidity of the environment, surface energy of the materials involved in the tribosystem, fatigue, and the presence of other particles may all play a smaller or bigger role in friction and wear.
Remark 2. Accurate and repeatable friction coefficient measurement is a challenging task, not only due to the non-equilibrium stochastic nature of the
sliding friction that leads to wear, but also due to the common lack of assessment (or at least presentation of assessment) of experimental uncertainty associated with the instrumentation. These uncertainties comprises several factors such as misalignment of transducers’ axis, force calibration, or voltage measurements, be- sides the measurand itself [78, 13, 61]. Furthermore, discussion on how common a steady-state of wear is can be found in [9].
2.5 Wear modelling
Wear modelling involves a few tasks, such as finding pressure distribution and deformation, depending on the intended approach. Initially, the contact problem needs to be described and solved, which may involve descriptions such as the Win- kler model [68], or by the use of half-space theory [3]. Next, a method to solve the equation must be employed, and again, a few options are available such as the Boundary Element Method, Fourier Transforms techniques, or the Finite Element Method. At last, the setting of the wear model is done, for which a common start- ing point is Archard’s wear equation. The chosen approach will define the amount of programming work; well-established methodologies with slight modifications are easily employed in practical problems by means of commercial software, while more state-of-the-art methods make use of more idealized problems to study meth- ods and assumptions themselves, demanding more or less computational and/or mathematical development by the researcher.
The fundamentals of contact modelling can be found in the well-known book Contact Mechanics by Johnson [43], and the numerical part can be found, for instance, in [93, 3]. A brief overview of the relevant concepts for this work is presented in the following subsections, followed by the literature review of scientific articles that inspired the approach detailed in the next section.
2.5.1 Contact modelling
A straightforward introduction to the formulation of contact is presented by [93]
and it is useful to transcribe part of it in this subsection. Consider a mass m under gravitation load and supported by a spring with stiffness k, as illustrated in Figure 4-a.
𝑢 ℎ
𝜆 𝜖
𝑎) 𝑏) 𝑐)
Figure 4: Point mass supported by spring system (a); illustration of the Lagrange Mul- tiplier Method (b); penalty spring due to penalty term (c). Adapted from [93].
The energy of the system, E, can be written as a function of the movement u as
E(u) = 1
2ku2 − mgu (4)
for which the extremum can be found by variation, i.e., applying the principle of minimum total potential energy, we have
δE(u) = kuδu − mgδu = 0. (5)
The second variation of E(u) results in k, which shows that the extremum of Eq. (4) is a minimum at u = mg/k. The restriction of the motion of the mass by the rigid support is expressed by
c(u) = h − u ≥ 0 (6)
which means that for c(u) > 0 there is no contact, and for c(u) = 0 the contact occurs. However, the variation δu is constrained by the contact surface, which from Eq. (6) yields that the virtual displacement is constrained and can only point in the upward direction, i.e., δu ≤ 0. Hence, using the variation in the variation form in Eq. (4) we have the following variational inequality
kuδu − mgδu ≥0 (7)
which characterize the solution of u. At this point, special methods can be em- ployed to solve the problem.
Recalling that once the contact has been established, a reaction force (R) will appear, the contact problem where the motion is restrained by Eq. (6) can be summarized as
c(u) ≥ 0, R ≤ 0 and Rc(u) = 0 (8)
which are known as Hertz–Signorini–Moreau conditions in contact mechanics. In the theory of optimization, such conditions coincide with Kuhn–Tucker comple- mentary conditions. The development above can be performed to include friction in the contact as well. In classical contact mechanics, the reaction force is assumed to be negative, i.e., only compression is considered.
Lagrange Multiplier Method. The contact problem constrained by an inequal- ity can be approached by using the method of Lagrange multipliers (Figure 4-b).
In this method, a term containing the constraint, which is equivalent to a reaction force, is added to the energy equation
δE(u, λ) = kuδu − mgδu + λc(u) = 0. (9) The variation of equation Eq. (9) results in two equations, since δu and δλ can be varied independently, one representing the equilibrium and the other rep- resenting the kinematic constraint. Thus, the variation is no longer restricted and the Lagrange multiplier λ can be solved, since it is equivalent to the reaction force λ = kh − mg = R. The solution of λ needs to fulfill the requirement of R ≤ 0, otherwise the contact is inactive and the solution is simply computed from Eq. (5).
Penalty Method. In this approach, a so called penalty term is added to the energy equation for an active constraint
δE(u, λ) = kuδu − mgδu + λc(u) = 0. (10) The penalty parameter can be interpreted as a spring stiffness in the contact interface. The variation of the equation results in
kuδu − mgδu − c(u)δu = 0 (11)
from which the solution and hence, the constraint are
u= mg+ h
k+ c(u) = h − u = kh − mg
k+ . (12)
Recalling that for contact we have mg > kh, this implies that penetration (depending on the penalty parameter) of the point mass into the rigid support occurs, which is physically equivalent to the compression of the spring (Figure4- c). Mathematically, the correct solution is approached as → ∞, which implies c(u) → 0, hence fulfilling the constraint equation.
It is important to highlight that both formulations exist in the normal and in the tangential direction of contact surface, which allows their computation in cases of rough contact with friction.
2.5.2 Finite Element Method
Partial differential equations (PDEs) are mathematical representations of the laws of physics for space- and time-dependent situations. Analytic solutions are avail- able only for a handful of cases, so approximations of the equations are performed for general cases.
Figure 5: A plate under uniaxial tension with a hole; finite element (FE) mesh of the model (middle) with finer (smaller) elements in regions where the fields and errors are expected to be high; typical visualization of a simulation. Adapted from [25].
The Finite Element Method (FEM) is a numerical method to compute approx- imations of PDEs; it generates solutions of boundary value problems in a weak
form, i.e., integral expressions obtained by the minimization of energy functions, using variational principles, or residual methods.
The approximation is based on the division (discretization) of the continuum domain into geometric entities named elements containing vertices named nodes, forming a mesh. The unknown field of interest - for instance, displacement field in a structural analysis - is piece-wise approximated by means of shape functions defined over each element and expressed in terms of nodal values. In other words, the field is built in the interior of elements by interpolation of the nodal values, which are obtained by solving the discretized problem. Figure 5 illustrates how a continuum domain is discretized and the result one can expect to see after a simulation, making use of a symmetry of the particular case to reduce the com- putational costs. In this manner the FEM allows to solve complex engineering problems in several areas of interest, which includes objects in arbitrary shape.
A general contact algorithm in a FE software consists of an initial search algorithm that identifies the penetration and then the satisfaction of the kinematic conditions, as described in the previous subsection. Commercial software offer such methods.
2.5.3 Archard’s wear equation
Once the contact problem is solved and information on contact pressure distribu- tion is available, the variable can be used as input to wear models. One of the most famous mathematical wear description of a sliding system was developed by Holm and Archard in the 1950s and is commonly known as Archard’s wear equation (or Archard’s law). It is arguably the most fundamental yet empirical model of wear and has become a standard in the field of tribology.
𝑣 𝑣
2𝑎
Figure 6: Wear occurs when a fragment of material detaches from an asperity during contact. The shape is approximated as a hemisphere with radius a. Adapted from [86].
Considering adhesive type of wear, the model assumes that sliding spherical asperities deform plastically in the contact, which have yielding limit as a hardness H. Hence, the mean contact normal load can be written as
dW = Hπa2 (13)
where a is the radius of the contact area. It is assumed that a hemisphere-shaped particle (Figure 6) with volume 2πa3/3 gets loose after sliding a distance of 2a with a probability of P . The wear volume per sliding distance is then
dV
ds = P2πa3
6a = P πa2
3 = P dW
3H . (14)
By integrating with respect to the sliding distance, the total wear volume V , can be expressed as function of the load, sliding distance and hardness as
V = KW s
H . (15)
Notice that we introduced K = P/3, which is known as the dimensionless wear coefficient, commonly obtained experimentally and used to characterize the material wear [91]. An analogous development can be made to obtain a similar expression to Eq. (15) for abrasive wear, where the wear volume is analyzed from a conical penetration approach. K can be interpreted in many ways: as the probability that the encounter of two asperities will originate a debris - and so, notice that Archard’s model is essentially probabilistic -, as the ratio of the volume worn to the volume deformed, as a factor reflecting the inefficiencies associated with the various processes involved in generating wear particles and others [74, 35].
Introducing the dimensional wear coefficient k = K/H (which is named specific wear rate[92]) and dividing Equation (15) by the area, an expression to calculate wear locally if the contact pressure is known is obtained:
h= kps (16)
where h is wear depth at the surface, p is the contact pressure, and s is the sliding distance. Furthermore, wear is a dynamic process and a wear rate equation can be simply obtained by deriving Eq. (16) with respect to time, getting the slip velocity vs in the equation. In simple terms, if the contact pressure and velocity vary in time, we write
dh(t)
dt = kp(t)vs(t). (17)
Evidently, Eq. (17) is a simplification of a process that was discussed (sub- section 2.4) to be a function of many factors. Nevertheless, it is a continuous description of the process that can be very useful as a starting point in numerical modelling. Clearly, Eq. (17) can be made more complex by adding effects of other phenomena such as presence of lubricants or oxidation, for instance.
2.5.4 Wear modelling in the literature
In the reviewed literature, Johansson [42] was one of the first authors to eval- uate wear in a numerical simulation; using the Finite Element Method (FEM), the contact pressure between two elastic bodies as the material is worn out by fretting4 was studied, presenting the governing equations in a discretized version and calculating the wear locally using Archard’s law. Debris or third body are not accounted for, i.e., the particles disappear when formed from the geometry.
As mentioned in the beginning of this section, the implementation of the wear model can be performed in a variety of ways, requiring more or less mathematical development by the researcher; a few examples are: Sawyer [77] presents closed
4Process in which two contacting bodies are subjected to slight relative vibratory motion, e.g.
riveted joints.
form expressions for contact loads and wear for a circular cam on flat faces fol- lower; Cantizano et al. [14] modelled different mechanisms of wear by means of definition of micromechanical law for contact and wear; in the work of Molinari et al. [57], Archard’s law was generalized to some extent to include features such as hardness dependence on temperature, allowing to capture wear regime tran- sitions; P´aczel and Mr´oz [63] discuss a variational formulation of contact shape evolution associated with wear. A few works have treated the modelling through what could be classified as a mechanistic approach: ratchetting wear mechanism for pin-on-disc is computed by Yan et al. [94]; sub-surface cyclic ratchetting plas- tic strain accumulations is the main idea in the wear model of Boher et al. [11].
Approaches to model fatigue in terms of accumulation of shear strain and layers is presented by Franklin et al. [24], and wear through crack growth involving linear elastic fracture mechanics and contact theory can be found, for instance, in the work of Ko et al. [46]. Ireman et al. [40] provide mathematical formulation of FE algorithms for two thermoelastic bodies in frictional wearing contact, presenting the governing equations in strong and weak form. Consideration of debris have been simulated in a semi analytic approach for fretting wear by Done et al. [20].
Different methods can be found for the definition and computation of the contact problem, often in an attempt to reduce computational costs: P˜odra and Andersson [68] use the Winkler model to evaluate simple cases; Flodin and Ander- sson [23] used it for wear simulation of helical gears, modifying Archard’s equation using Hertzian solution; Peigney [64] develops a minimization approach to obtain the worn-out geometry. Sharif et al. [80] modelled wear of worm gear under elasto- hydrodynamic lubrication, using finite differences method and adapting Archard’s equation to include lubricant film thickness. Sfantos and Aliabadi [79] used the Boundary Element Method to model contact and simulate wear. Andersson et al.
[4] presents an implementation using Direct Convolution and Fast Fourier Trans- form to approach the contact problem and Archard’s wear equation to compute time dependent wear in a sphere on flat contact. Asymptoptic modelling is per- formed by Argatov and Tato [5]. Recently, Cavalieri et al. [15] used a specialized FEM with a so called mortar algorithm to study wear in internal combustion engine valves.
Early works in wear modelling often would not account for geometry change, using contact pressures and slip from an initial wear cycle, leading to unreliable results. The coupled evolution of contact and geometry, however, is of uttermost importance; most approaches on wear modelling, since the development of com- putational algorithms, became more and more accessible and have been based on evolving contact conditions in different cases: the contact geometry varies gradu- ally in an iterative procedure in which the relevant results - such as pressure and time/sliding interval - are progressively computed at each iteration and used in a wear model, which defines the change in the geometry.
FEM & Commercial software. The general idea in modelling through evo- lution of contact conditions has been applied in may fields and problems. P˜odra and Andersson [67] used FE analysis with commercial software ANSYS (i.e., a general-purpose FE software) to obtain the contact solution and simulate sliding wear by means of Archard’s model of the wear depth with respect to sliding dis- tance, of a hemispherically tipped pin sliding against a plane. The basic algorithm consists of solving a series of structural static analysis with the geometry being
updated due to the calculated wear depth. Issues regarding the size of wear incre- ment are mentioned. The authors also presented a version which included surface topography on the wear model by means of linearisation of the Abbot curve [66].
¨Oqvist [62] performed a similar study, but also accounting for elastic deformation in the contact zone.
Yan et al. [94] present a 2D computational approach to predict the sliding wear in terms of accumulated ratcheting strain in a loaded spherical pin con- tacting a rotating disc (pin-on-disc), idealizing the problem as a 2D plane strain approximation, thus reducing computational costs. By using a static 3D simula- tion, a conversion factor for the wear rate associated with a 3D sliding contact is introduced.
Hegadekatte et al. [35] developed a wear processor along with commercial FE package ABAQUS for dry sliding wear, using Euler integration for the wear over the sliding distance, so the simulation also works as a series of static simulations with updated geometry; the author highlights care taken with the wear direction:
nodes are shifted inward a surface normal direction, defined by the edges of the adjacent finite elements. Later, Hegadekatte et al. [34] used their algorithm to simulate the wear of a micro-planetary silicon gear train.
McColl et al. [55] employ a 2D FE model in ABAQUS to study frictional contact of a cylinder-on-flat fretting test, using a modified version of Archard’s equation for sliding wear. The FE mesh at the contact is as fine as 10 micrometers.
Gonzales et al. [29] simulated the experimental pin on disc tests of an Al–Li/SiC composite using ABAQUS, using an isotropic thermo-elastoplastic, but a simpli- fied 2D analysis.
Figure 7: Geometric model of disc brake assembly on which wear of the pad was modelled and numerically studied; Surface profile of brake pad for wear simulation. From [84].
S¨oderberg and Andersson [84] studied the pad-to-rotor interface of car disc brakes as a conformal dry sliding contact (Figure 7) to calculate pressure and wear; they used ANSYS and the macro scale, phenomenological Archard’s wear law. Only wear at the pad was considered and the sliding distance was calculated from a rigid displacement. The pad is assumed to wear perpendicular to the mating rotor surface (which is rigid); this is an important and interesting assumption, because it facilitates the geometry update algorithm; they do make a remark on this matter: “A more correct procedure would be to assume that the pad surface is worn in the inwards direction of its local surface normal, which may vary over time due to uneven wear over the contact surface. The simplification is justified by the fact that the two contact surfaces are initially flat and conformal, and that the
expected wear depth is much smaller than the dimensions of the contact surface.”.
They use the same approach as P˜odra and Andersson [67], which is similar to Hegadekatte et al. [35] and ¨Oqvist [62], i.e., it works by looping through a series of static structural analyses of the FE model, each with an updated geometry of the pad surface, but without FE re-meshing between steps. Another factor that may influence the accuracy of the result is the step length during the numerical integration of the wear equation. The accuracy of the solution obtained can be evaluated by comparing solutions with different step lengths.
Numerous simulations. In general, wear occurs slowly, but the linearisation of the process using Archard’s wear equation and the approach of incrementally simulating the coupled contact and geometry can lead to high computational costs when performing numerous simulations. Therefore, it becomes desirable to per- form extrapolations on the simulations, raising issues on the reliability of the results, briefly highlighted by Dickrell et al. [19]. The approach on this matter done by Kim et al. [45] in simulating a block-on-ring experiment is done by assum- ing that the contact pressure of a single FE analysis is constant for several cycles:
initially, a single oscillation cycle is discretized into Ns number of steps, on which the wear model is computed; then, it is assumed that No oscillation cycles behave in the same way. Next, after performing NF E FE analyses, the results would be representing NF E× No experiments. Finally, assuming that at this stage a steady progression of wear occurs, the wear depths are linearly extrapolated Nc times, corresponding to NF E × No× Nc cycles.
Mukras et al. [58] uses an extrapolation scheme by including an “extrapola- tion factor” to the incremental wear depth calculation, projecting the wear of a cycle to several cycles. A condition is placed (and justified) such that the overall smoothness of the pressure distribution should remain. The approach used by the authors is discussed again in the next section. The authors also proposed a procedure to use an adaptive factor in order to try to find a optimum value for each cycle, in terms of deformation and contact characteristics.
Chongyi et al. [16] studied a wheel/rail wear profile using Archard’s law in a spinned 2D axisymmetric model in ABAQUS software, using limitations or smoothing in order to avoid unrealistic scenarios, i.e., “to get better approxima- tions of the continuous wear process with a discrete sequence of profile updates”.
In order to do that, a cubic spline interpolation algorithm is applied on the wear distribution before starting a new iteration of the process. Another interesting aspect is the use of bilinear elastic–plastic constitutive material model.
Other interesting works. Rezaei et al. [73] studied polymeric composite journal bearings with ABAQUS, modelling an orthotropic material and validat- ing through experiments. In a subsequent work, Rezaei et al. [72] simulated wear following the general idea of evolving contact conditions, which they named “adap- tive wear simulation”. The wear direction, i.e., sweeping of FE nodes in the inward normal direction was detailed following the same idea as Hegadekatte et al. [35].
Rezaei et al. [72] was presumably the first to simulate wear using a non-isotropic material model. They also applied a pressure and velocity dependent wear coeffi- cient. In a flat on flat sliding simulation they applied a constant pressure on the contacting surface, instead of taking it as output.
K´onya and V´aradi [47] studied wear of polymer-steel sliding pin on disc pair, considering temperature and time-dependent behavior of polymer materials; in
their algorithm, the thermal expansion of the polymeric pin is considered in the wear increment, having an opposite effect to the wear until a steady state is reached. Mart´ınez et al. [53,52] also performed wear modelling in polymer-metal contact, using an experiment-fitting power law for the wear model, performing different simulations to extrapolate the model for large cycles. However, linear elastic material model was used.
When using general-purpose FE software, the possibility of using coupled thermo–mechanical analysis is eased, which can be used for instance to include the effects of frictional heating. Shen et al. [82] performed a 2D thermo-mechanical simulation in ABAQUS of a spherical plain bearing with a self lubricating compos- ite liner. The consideration of the thermal analysis was performed by translating the temperature distribution into body loads. Lengiewicz and Stupkiewicz [50] de- veloped a general periodic pin-on-flat wear problem, considering wear on both con- tacting bodies. Gui et al. [32] studied thermo-mechanical simulations in a brake system. Zhang et al. [95] studied the tribological properties of self-lubricated polymer steel laminated composites and simulated the thermal expansion of the solid lubricants. In their experiments, they calculated the wear coefficient bas- ing on wear marks and geometrical features. In their simulations, the transient temperature field was studied with the account of a few assumptions so that the conservative law of energy was applicable; material properties were independent of temperature.
The general methodology with slight modifications has been applied to many different problems: cam-followers works [37,38,85], wear of artificial cervical discs [8], wear in total disc replacements [71], wear and tribofilm growth [3], machining tools [70], and press hardening tools [18].
3 Research Gap
In section 2, it was shown that there are many events occurring in a tribological contact which are dependent on many physical phenomena. Numerical modelling is a powerful tool that can be used to try to model such events and make pre- dictions possible. It was seen in the literature review that many simplifications must be performed to make simulations feasible to be analyzed; depending on the approach, published studies focus on a practical problem using more established methodologies or evaluate the numerical method itself.
In context of dry sliding, frictional contact and wear, the main goal of this thesis is to analyze a practical problem. The objective can be stated as: the devel- opment of a numerical model and method to predict contact pressure distribution, wear and contact geometry evolution of a linearly reciprocating block in contact with a flat plane. In the reviewed literature, no work on modelling was found to approach this particular problem, although the general methodology of progressive computation of contact parameters in a iterative fashion is applicable to it.
Contact pressure and local wear during sliding is not something easily mon- itored during an experiment, thus, real time information about these entities’
distributions in such cases are a valuable result intended from this work. Such information would result in a better understanding of the linearly reciprocating dry sliding contact. The scope of the thesis is subjected to the time-span available, thus defined to be limited to simulating the above mentioned entities considering only solid mechanics and contact mechanics physics, intended to be a first stage model of the proposed case.
The questions expected to be answered regard not only how the evolution of the entities occur, but also how to set up a methodology using available software. The development of the model is described in section4along with the assumptions and approximations, which can be critically analyzed with the theoretical knowledge presented in the next section.
4 Model & Methods
In this section, the tools used and the details on the development of the model are presented. Initially, the tools and model problem are presented along with the relevant aspects and considerations used. Secondly, the employed methodology for wear and geometry update inspired by the literature review is presented.
4.1 Software
Each approach in contact mechanics and wear modelling may have its advantages and disadvantages in terms of accuracy and computational cost [58]. A common approach is to use the Finite Element Method, by means of general-purpose com- mercial software. Although it has the disadvantage of being time-consuming, FEM is a powerful tool to solve contact problems and it is the approach used in this work. It is important to highlight that such approach involves explicitly modelling of the geometry, setting up the materials, physics, FE mesh and so on to obtain the solution at the contact.
In this work, COMSOL Multiphysics® 5.3 and LiveLink™ for MATLAB® are the software used. COMSOL Multiphysics® is a finite element analysis simulation software, which allows coupling of systems and physics [17]. Among the several modules, the Structural Mechanics Module allows the static and dynamic analysis of mechanical systems, i.e., a wide range of analysis types, including time studies in stationary or transient regimes, eigenmode/modal analysis, parametric, quasi- static, frequency-response, buckling, and prestressed. LiveLink for MATLAB ex- tends the COMSOL modeling environment with an interface between COMSOL Multiphysics® and MATLAB®, which may be used to control the software ex- ternally, i.e., through the scripting programming in the MATLAB environment, handling the COMSOL model in preprocessing, manipulation, and postprocessing.
The simulation workflow in COMSOL basically goes through the following steps: set up of the model environment, creation of geometric objects, specification of material properties, definition of physics and boundary conditions, creation of FE mesh, computation, and post-processing. Evidently, the settings of each step shall be performed in order to best model the problem in hand, which will be discussed in the next subsection.
4.2 Model problem
3D modelling has the advantage of considering a correct geometry (dimensions, surface topography, degrees of freedom) and ability to calculate stresses and strains in the entire bodies. However, 3D-FE models demand very high computational costs and processing times, which turns out to be a disadvantage when method- ologies are being studied and parameters are being varied.
Considering the problem is a first stage model and the methodology for wear is to be studied, which involves experimenting with different numerical parameters, the block in linear reciprocating motion in contact with a flat plane is approxi- mated by the use of a 2D space dimension. Among the interfaces available in the Structural Mechanics Module, the Solid Mechanics interface is used, which can provide information on displacements, stresses, and strains of the simulation.
The 2D approximation can be performed using plane strain or plane stress conditions. In plane stress, the normal stress and the shear stresses directed perpendicular to the plane of study are considered to be zero. Hence, this condition is appropriate for bodies whose thickness is small. In plane strain, the normal strain and the shear strains on the plane of study are assumed to be zero, which is appropriate when the out of plane dimension is much larger than the other two [27]. In the modelling of this work, the plane strain condition is used. Additionally, no inertial terms are included in the time domain analysis.
The linearly reciprocating block in contact to a flat plane is modelled, essen- tially, as two bodies (the block and the counter-surface, Figure 8) and simulated in two study steps in COMSOL: a Stationary and a Time Dependent. Solving the problem in these two steps is - besides a good approach in assuring the com- putation will likely be successful at the end (in terms of solver convergence) - in accordance to the actual experimental case; initially, the sample block is positioned and the static load is slowly applied before the movement starts.
From the “point of view” of the software, every surface of an object can come in contact to every other surface, so defining a Contact Pair under Definitions is important for better computational efficiency, because then the software can limit its analysis operation of contact search. Contact is typically seen as one-sided rigid constraint; the contact algorithm in COMSOL constrains the destination boundaries so that they do not penetrate the source boundaries. As a general rule, one should make the stiffer object the source object or, if the stiffnesses are similar, make the concave side the source.
The contact feature5 makes the studies geometrically non-linear6. Commonly, strain calculations are treated in a geometrical linear fashion, i.e., the equilibrium is formulated in the underformed configuration and is not updated with deforma- tion. This implies small errors, but which are not so big to the point of justifying the implementation of a more mathematically complex equilibrium model, in gen- eral cases. In the case of contact physics, however, it cannot be used in a linear context [87]. Fortunately, COMSOL is equipped to deal with such cases.
Table 1: Material properties in the model, Thordon’s properties obtained from [90] and steel data from COMSOL’s material library.
Block Counter-surface Young’s modulus [GPa] 2.41 205
Poisson’s ratio [-] 0.39 0.28
Density [kg/m3] 1400 7850
Once the geometry is defined, the materials can be added to the components.
The Thordon Bearings’ product ThorPlas®, a homogeneous, self-lubricating ther- moplastic material designed to operate in high pressure environments such as hydro-turbine wicket gate in hydropower plants [89], is the material used for the block. The reason for this choice is due to internal interest on this material at the university. The material chosen for the counter-surface is the steel AISI 4340.
5Applicable only to solid elements in COMSOL.
6When geometric nonlinearity is included in structural analysis, the COMSOL automatically makes a distinction between Material and Spatial frames. The material frame corresponds to the undeformed configuration, while the spatial frame corresponds to the deformed configuration.
The important material properties are shown in Table 1. For simplicity a linear elastic model was used as material model for the block.
With the aid of Figure 8, the numbered aspects and boundary conditions are detailed for the contact analysis; the bottom region of the counter-surface is sup- ported by the machine. Hence, a boundary condition of fixed constraint 1 is applied to the bottom of the counter-surface, which means that all translations and rotations are zero on that boundary.
y x
1
2 5
4
3
Figure 8: Left: illustration of block on plate experiment (from [28]); right: modelling of experiment - refer to text for numbered explanation.
The counter-surface geometry is built in such a way that allows an appropriate preparation for the simulation in terms of FE mesh. A division of the domain is performed 2 for meshing purposes, although it consists of the same body.
Recalling that the FEM approximates the solution by means of discretization of the domain of interest, the accuracy of the solution is linked to the mesh size, and so, it is important to have a sufficiently fine discretized mesh of the domain7. Figure9 shows the FE mesh for the counter-surface.
Figure 9: FE mesh for the counter surface showing finer refinement towards region where contact with the block occurs.
The modelling of the contact was performed in such a way that the bodies are not in contact in the initial configuration, i.e., a gap 3 exists between the block and the counter-surface. The reason is that when wear occurs, regions previously in contact will not be in contact anymore, so there will be a creation of a gap. Thus, the modelling needs to be able to deal with a contact modelling of bodies which might not be in contact at the start of the simulation. As a consequence, the block is initially unconstrained and there will be possible rigid body displacements, which often causes the solver to fail. This issue is avoided in the load and displacement boundary conditions.
7Theoretically, as the mesh size approaches zero, the solution approaches the exact solution, but limitations with regards to finite computational resources and time exist. Thus, the goal of a simulation is to minimize the difference (“error”) between the exact and the approximated solution, according to some accepted tolerance criterion.
A rigid connector feature is used on the top of the block 4 . Rigid connectors are boundary conditions for modelling rigid regions; the selected boundaries will behave as a single rigid object. In reality, the load on top of the block is applied by the steel parts of the machine, and thus the rigid connector is appropriate in this case not only because the stiffness of the steel is much larger than the block’s material, but also and more importantly because it simplifies the simulation. Next, a functionality of an applied (parameterized) force is added to the rigid connector on the top of the block.
Another rigid connector 5 is used in the upper parts of the block’s side boundaries - steel parts of the machine contact the sample block at those locations in the actual case - in order to apply the displacement for the Time-dependent study. The displacement was imposed using an expression of a wave function globally defined to represent the oscillatory linear motion, which can be expressed as a function of time according to
u(t) = A sin (ωt + φ) (18)
where A is the amplitude of the displacement - half the stroke - ω is the angular frequency of the motion, i.e., in [rad/s], and φ is the phase, which is useful to define the starting position of the block in the simulation (providing initial values of displacement field are also defined).
In a two step simulation, the initial values of the second step should be con- sistent with the final solution of the first step. Hence, a global parameter t = 0 is defined to be used in the Stationary study step, which results in no displacements from the displacement function dependent on t.
Another important aspect refers to the initially unconstrained problem; a Spring Foundation feature is added to the boundaries of both 4 and 5 , which are in essence, a temporary weak spring during the beginning of the simulation.
With the springs present, the rigid body motions are constrained, hence the prob- lem can be solved. In order for them to not affect the actual configuration of the problem, a parameter P ar ranging from 0 to 1 is used to define a stabilizing spring’s stiffness ku as
ku = k0(1 − P ar) × 2−(P ar×10) (19) which results in the stiffness reducing to zero as the parameter reaches 1. k0
must be chosen in such a way that the spring force balances the external load at a sufficiently small displacement; “a too weak spring will give a too large initial overclosure of the contact boundaries, and a too stiff spring might influence the solution too much” [17]. In this modelling, the value was chosen at 106N/m. This approach avoids the rigid body displacement that could occur on the block.
The load applied at 4 is also parameterized. In order to “help” the soft- ware converge better, it is useful to start the iteration close enough to the solu- tion (Newton-Raphson or similar methods are used); this is done by ramping the boundary condition, in this case, the compressive load. So, a small load can be applied and once the simulation has converged, we use the solution as a starting point when iterating for a higher load until the desired load has been reached.
This is done with the parameterization Par × Load and in adding a auxiliary sweep in the solver of the software. Thus, the same parameter Par is used for
