Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 50, 2, pp. 509-530, Warsaw 2012 50th Anniversary of JTAM THE MODELLING OF LOADING, UNLOADING AND RELOADING OF THE ELASTIC-PLASTIC CONTACT OF ROUGH SURFACES Robert Kostek University of Technology and Life Science in Bydgoszcz, Faculty of Mechanical Engineering, Bydgoszcz, Poland; e-mail: robertkostek@o2.pl This paper presents a non-linearmathematical model describing the hy- steresis of dry contact of rough surfaces loaded in thenormaldirection.A number of characteristic phenomena appear during the initial loading, unloading and reloading of the contact interface. The contact is most flexible during the first loading, and its plastic deflection is significant, whereas the characteristic curve of unloading is stiffer and ismainly ela- stic in nature. The reloading characteristic is similar to the unloading characteristic, and thus a hysteresis loop is observed. The non-linear model of contact presented in this paper, which has been validated by experimental data, describes the aforementionedphenomena.Themodel is dedicated to the finite element method (FEM), the discrete element method (DEM), and simulations of multi-body systems (MBS), where accurate and effective contactmodels are preferred. The presentation of themodel is supported by examples of applications to contact problems. Key words: contact model, elastic-plastic, rough surface, hysteresis 1. Introduction Contact interfaces often exist in machines (e.g., sideways). The contact in- terfaces operate under dynamic loading and are thus repeatedly loaded and unloaded. Contact deflections are particularly important in the case of preci- sionmachines, such as machine tools, measuring equipment andmicroscopes. The specific physical phenomena that occur at the contact of two rough surfaces are strongly non-linear. The maximum normal force influences the residual deflection of the contact (Etsion et al., 2005; Kadin et al., 2006). The plastic deflection of a virgin contact is significant (Goerke and Willner, 2008; 510 R. Kostek Tworzydlo et al., 1998; Vu-Quoc and Zhang, 1999), whereas a stabilised con- tact is stiffer, and its deflection ismostly elastic in nature (Skrodzewicz, 2003). Moreover, the hysteresis loops of stabilised contact joints are not very sensiti- ve to the frequency of external forces, within the frequency range 0.1-2400Hz (Gutowski and Skrodzewicz, 2001; Kostek, 2005). The specific phenomena make it difficult to model the contact; thus, contact models are dedicated to specific applications. In tribology, accurate and complex models are pre- ferred. They are based on statistical modelling (Kadin et al., 2006), fractal theory (Goerke and Willner, 2008; Pei et al., 2005), FEM (Tworzydlo et al., 1998; Pei etal, 2005), and modelling dislocations within crystals (Nicola etal, 2007). Analysis of thesemodels provides valuable information, but simulations using suchmodels are computationally expensive. In contrast, simple and fast models based on linear (Kruggel-Emden etal, 2008; Munjiza, 2004) and non- linear formulae (Grudziński andKostek, 2007;Hess andSoom, 1991;Hunt and Crossley, 1975; Kostek, 2005; Kruggel-Emden etal, 2008; Martins etal, 1990; Michalczyk, 2008;PomboandAmbrósio, 2008; Schiehlen etal, 2006) areprefer- red in computational mechanics and dynamics. Thesemodels are usually used to describe collisions between bodies, vibrations, pressure distributions and restitution coefficients, which are coupled with the contact phenomena (see, e.g., Lifshitz and Kolsky, 1963; Oden andMartins, 1985; Wriggers, 2006). The contact interfaces of rough surfaces have been investigated by resear- chers, and contact models describing loading and unloading seem to be the most valuable. Tworzydlo et al. (1998) presented asperity-based constitutive models of contact interfaces. Skrodzewicz (2003) presented non-linear mathe- maticalmodels describinghysteresis loops of stabilised contacts loaded byhar- monic forces.Kim et al. (2004)worked onultrasonic experiments characterised by the elasto-plastic planar contact. Jones (2004) presented a model of a ro- ugh surfacewith columnar and spherical asperities. Pei et al. (2005) presented an FEM study of the contact between a rigid plane and an elasto-plastic solid with a self-affine fractal surface.Kadin et al. (2006) presented a statisticalmo- del of the unloading of the elasto-plastic contact of rough surfaces.Goerke and Willner (2008) worked on a numerical model that describes the elasto-plastic normal contact of isotropic fractal surfaces. Kostek and Konowalski (2008) presented a non-linear model describing the loading and unloading contact of two rough surfaces. Several articles have addressed elasto-plastic Hertzian contact in the Greenwood-Williamsonmodel (Archard, 1957; Etsion et al., 2005; Greenwood and Williamson, 1966; Jamari and Schipper, 2007; Li and Gu, 2009; Pullen andWilliamson, 1972; Schiehlen et al., 2006; Vu-Quoc and Zhang, 1999;Woo The modelling of loading, unloading and reloading... 511 and Thomas, 1980), which can be relevant to the contact of asperities. Thus, these articles have beenmentioned in this paper. As is evident from the literature search, there aremanymodels of contact. However, a model that has been validated by experimental data and that is able to precisely describe the initial loading, unloading and reloading in the normal direction of a planar contact of rough surfaces has not been found. Thus, there is a need to present a suitable model for that case. 2. Formulation of the mathematical model The interface of two bodies in a planar contact, which is shown in Fig.1a, is modelled in this Section.The interface can bemodelledwith a large number of springs and plastic bodies (Saint-Venant elements) (Fig.1b), which represent interacting roughness asperities. It is assumed that the mechanical properties of the contact interface are macroscopically identical over the entire contact area.Thus, after homogenisation (Fig.1c), the contact zone canbemodelledas Fig. 1. Schematic of the contact interface (a) and its physical models (b), (c) a homogenous area that possesses specific elasto-plastic properties. Deflection of the contact zone can be expressed as a function of loading, where the input signal is the normal contact pressure p, and the output signal is the contact deflection δ. The experimental and analytical investigations presented previo- usly (Mikic, 1971, 1974; Tworzydlo et al., 1998) show that the interactions within a dry contact exhibit elastic, plastic and elasto-plastic nature (Fig.2d), and thus the contact model consists of three elements (Fig.2a). The elastic normal deflection is described by equation (2.1)1 (Martins et al., 1990), and the plastic normal deflection is expressed by similar equation (2,1)2 as follows δe = eep me δp = epp mp max (2.1) where δe denotes the normal elastic deflection of the interface, δp represents the normal plastic deflection of the interface, p is the normal contact pressure, 512 R. Kostek pmax is the maximum normal contact pressure for a specified period of time pmax =max(p0,p1,p2, . . . ,pn), and ee, ep, me and mp are the parameters of the interface. The non-linear elastic element models the reversible part of the contact deflection, whereas the non-linear plastic element models the irrever- sible part of the contact deflection. Themost important feature of the elastic element is that its deflection does not depend on the previous loading history, only on the present contact pressure p, whereas the plastic element deflection is determined by themaximummagnitude of the loading pressure pmax. The plastic contact deflection increases only if themaximal contact pressure pmax increases. The two elements are able to describe pure elastic and pure plastic deflections (Kostek and Konowalski, 2008; Nardin et al., 2003). Fig. 2. A load-deflection characteristic (d) of the contactmodel (a) obtained for the presented input (c) and output signals (b) and validated against experimental results (see Tworzydlo et al., 1998; Fig.20a) The application of the elasto-plastic element provides an opportunity to model thehysteresis of apreloadedcontact andthedissipationof energywithin this type of contact. This element can bemodelled as a parallel connection of the elastic and plastic bodies (Fig. 3 a), described by equations (2.2). The pla- The modelling of loading, unloading and reloading... 513 Fig. 3. A load-deflection characteristic (d) of the elasto-plastic element (a) obtained for the presented time histories (b), (c) stic bodymodels the dissipation of energy and damping, and thus its pressure is related to the sign of its velocity of deflection. It should bementioned here that the damping pressure is insensitive to themagnitude of the speed, which reflects the plastic nature of dissipation energy. In general, the elasto-plastic element is described by the following equations pep−e = ( δep eep−e )1/mep−e pep−p = sgn (dδep dt )( δep eep−p )1/mep−p (2.2) and p(t)= pep−e+pep−p = ( δep eep−e )1/mep−e + sgn (dδep dt )( δep eep−p )1/mep−p (2.3) where δep denotes the normal elasto-plastic deflection of the interface, pep−e is the normal contact pressure in the elastic element, pep−p is the normal contact pressure in the plastic element, t is time, and eep−e, eep−p,mep−e and mep−p are interface parameters. In this case, the input signals are δep and dδep/dt, 514 R. Kostek whereas the output signal is the normal contact pressure p. First, during the initial loading from point A to point B, both pep−e and pep−p are positive (dδep/dt > 0), and thus these two elements act against the normal contact pressure p. In consequence, the elasto-plastic element follows the loading cu- rve (Fig.3d). Next, the unloading from point B to C does not change the deflection of the elastic element (dδep/dt = 0) (Fig.3b,d), but only changes the normal contact pressure in the plastic element pep−p, which is because the plastic element stops the changes of deflection, which is typical for sys- tems with dry friction. The contact pressure during this period drops from pB to pC, due to a drop of the external force. The contact pressure being in the range from pC to pB is not able to change the deflection of the elastic element. Then, the unloading from point C to point D reduces the elasto- plastic deflection of the interface δep. At this period, the plastic element acts against the elastic element (dδep/dt < 0), and the elasto-plastic element fol- lows the unloading curve. Thus, some hysteresis loop is observed. Finally, the reloading from point D to point F follows the loading curve (dδep/dt > 0). This description clearly shows that the elasto-plastic element has three possi- ble stages: dδep/dt> 0, dδep/dt=0, and dδep/dt< 0. If mep−e and mep−p are equal, which makes the formulae simpler, then formula (2.3) can be written for the three stages as follows if dδep dt > 0 then δep = eep−lp mep if ( δep eep−e ) 1 m ep−e + ( δep eep−p ) 1 m ep−p >p> ( δep eep−e ) 1 m ep−e − ( δep eep−p ) 1 m ep−p then dδep dt =0 if dδep dt < 0 then δep = eep−unlp mep (2.4) where eep−l, eep−unl and mep denote the interface parameters. The contact pressure can be treated as the input signal, and then it governs the switching between the three stages. An example for the stage dδep/dt = 0 presents equation (2.4)2. More information and pseudo-code are presented in the Ap- pendix. The elasto-plastic module can be interpreted as a Hertz-Stajerman model (Michalczyk, 2008). The contact model presented in Fig.2a consists of three modules: plastic, elasto-plastic and elastic. Deflection of themodel δ is the sum of three deflec- tions given by δ= δp+ δep+ δe (2.5) The modelling of loading, unloading and reloading... 515 where δ denotes the normal deflection of the contact. It is worth noting that themodel have only seven parameters: ee, ep,me,mp, eep−l, eep−unl,mep, and the contact behaviour described by the model is very similar to the experi- mental results (Fig.2d). The loading frompoint A to point B (Fig.2c) causes deflections of the three elements – plastic, elasto-plastic and elastic (Fig.2b) – which indicates that the contact model is flexible (Fig.2d). The unloading to point C changes the deflection in the elastic element only, whereas the unlo- ading below point C reduces the deflection of the elasto-plastic element. The elasto-plastic element follows the unloading curve,Eq. (2.4)3, at this time.Du- ring the reloading from point D to point E, the elastic and elasto-plastic ele- ments are deflected together, and the elasto-plastic element follows the loading curve Eq. (2.4)1. Exceeding the highest pressure so far, pmax−B at point E, makes the plastic deflection larger, and the contact becomesmore flexible aga- in. The three elements (modules) are then deflected together between points E and F. A similar sequence will be repeated during the next unloading- reloading cycle. The plastic nature of the dissipation is worth noting. Itmakes the model insensitive to the speed of deflection, because the damping force is insensitive to the magnitude of speed. This phenomenon is coupled with the insensibility of the dry contact hysteresis loops to the frequency of excitation (Gutowski and Skrodzewicz, 2001). Table 1.The contact data used for the simulations Parameters Magnitudes of contact Units interface parameters ep 0.534 µmMPa −mp mp 0.653 ee 1.148 µmMPa −me me 0.462 eep−l 0.302 µmMPa −mep eep−unl 0.523 µmMPa −mep mep 0.462 The parameters of themodel, which are presented inTable 1, were estima- ted from experimental results (Tworzydlo et al., 1998) using theMonte Carlo method. This method can be used to find the global minimum in a specific solution space. The contact pressure p was assumed to be an independent variable, whereas the contact deflection δ was assumed to be a dependent variable. The residual sum of squares was adopted as an objective function. Finally, thedifference between the computational results and the experimental 516 R. Kostek results was minimised, which leads to a good agreement between the compu- tational and experimental results (Fig.2d). Experimental investigations of the contact characteristics have beendescri- bed previously in literature, see, e.g., Cichowicz and Nowicki (1984), Goerke and Willner (2008), Grudziński and Jaroszewicz (2002), Gutowski and Skro- dzewicz (2001), Kostek (2005), Skrodzewicz (2003), Tworzydlo et al. (1998). The main difficulty associated with these investigations is the small contact deflection δ, which is a result of the small height of the roughness asperi- ties. To overcome this difficulty, special displacement transducers are applied and special experimental setups are designed, or common testingmachines are used. The experimental results used as a reference in this paper were obta- ined in Tworzydlo et al. (1998) for several disk specimens made of 1020 steel (Fig.4). The thickness of the specimens was 3.175mm and the diameter was 25.4mm. The contacting surfaces of the specimens were glass-bead blasted, which introduced roughness. The standard deviation of the profile heights was σ=4.91µm.A sandwich of the specimenswas first loaded, then unloaded and finally reloaded. Recording the time histories of the loading force and contact deflection provides the opportunity to determine the contact characteristics. A more detailed description of the experimental setup and procedure can be found in Tworzydlo et al. (1998) Fig. 4. Schematic of the experimental setup (Tworzydlo et al., 1998) 3. Application examples The model of contact presented in this paper, which was validated experi- mentally, has been applied to the modelling of interactions between rigid and elastic bodies. The parameters of the model are shown in Table 1. The first The modelling of loading, unloading and reloading... 517 example presents an interaction between rigid bodies that is typical for MBS systems.The second example presents the interaction between an elastic beam and a rigid flat. The last example presents the collision of elastics bodies that was simulated with the combined finite-discrete element method. 3.1. Rolling contact A Hertzian contact is a classic issue in contact mechanics and tribology that is coupled with the coefficient of rolling friction and the deflection of contact. The model used for this case consists of a rigid cylinder rolling on a rigid flat (Fig.5), and the elastic nature of the bodies is neglected (MBS). This model presents the influence of deflections of roughness asperities on the system. The contact zone has been modelled as a homogeneous layer, the constitutive relation of which has been described by themodel presented here (Fig.2). The contact pressure has been calculated for a number of points and then interpolated. In this way, the interaction between the bodies has been modelled. Fig. 5. Schematic of the rolling contact between a rigid cylinder and a rigid flat During rolling, the contact zone is first loaded and then unloaded; the two stages of the process are presented inFig.6. The loading process and the unlo- adingprocess followdifferentpaths: the A-B characteristic during loading and the B-C-D characteristic during unloading (Fig.2d), which introduces dissi- pation energy, rolling friction and residual deflection of contact δpl (Fig.6a). The dissipation energy is coupled with the contact hysteresis (different lo- ading and unloading characteristics). The amount of energy absorbed during loading is greater than the amount of energy released during unloading be- cause of plastic deformation asperities. In other words, the area under the 518 R. Kostek loading characteristic is greater than the area under the unloading characte- ristic, which reflects the dissipation energy. The energy loss is proportional to the area of the hysteresis loop A-B-C-D-A. Fig. 6. Characteristics of the rolling contact (a) and contact pressure distribution (b) obtained for the following data: D=100mm,L=100mm, F =611.7N The resulting distribution of contact pressure is asymmetrical (Fig.6b), which reflects contact hysteresis, dissipation energyanddeformations of rough- ness asperities (Stolarski, 2000, p.249). The asymmetry is a result of different loading and unloading characteristics. Thus, this model provides an opportu- nity to calculate the coefficient of rolling friction, which equals f =0.072mm. “It is quite obvious that resistance to the rolling of a wheel is greater on a rough surface than ona smooth one, but this aspect of the subjecthas received little analytical attention” (Stolarski, 2000, p.237). The second issue is the contact deflection. For the contact interface descri- bed here, the maximum contact deflection equals δmax = 6.07µm (Fig.6a), which is over ten timesmore than the contact deflection resulting from a clas- sical Hertzian contact of smooth bodies made of steel. As is evident from the results, the deflection of rough surfaces plays an important role in this system. 3.2. Prismatic beam on a rigid flat The second example involves a cube-shaped elastic steel beam resting on a rigid flat surface and loaded by a force F (Fig.7a). The elastic beamhas been modelled with twenty beam elements (Young’s modulus E = 2 · 105MPa), which is a simplification of themodel (Fig.7b). The contact zone between the beamand the rigid flat has beenmodelledwith the contactmodel presented in this paper (Table 1). The contact pressure p has been calculated for 21 nodes (Fig.7b), and then the contact pressure distribution between the nodes has been interpolated with a linear function. Consequently, the distributed load q has been interpolated with a linear function as well (Fig.7c). The loading The modelling of loading, unloading and reloading... 519 for a non-uniform space mesh can be calculated with this approach. Finally, the deflection of the beam due to contact pressure has been computed. The non-linearity of the contact presents some difficulties associated with the cal- culation of the contact pressure,whichwere solvedwith an iteration procedure and a gradual increase of the loading force F. Fig. 7. Schematic of the contact interface between an elastic beam (L=200mm, B=10mm,H =62.4mm) and a rigid flat (a), physical model of the interface (b), and interpolation of the load distribution between nodes (c) The contact interface has been loaded, unloaded and reloaded, and a hy- steresis of the interface has been obtained (Fig.8). During the initial loading, the maximum pressure and contact deflection are observed in the middle of the interface (Fig.9). As the loading force F increases and the beam distorts, the beam ends are gradually unloaded. The unloading process is different from the loading process. The contact pressure in the middle part of the beam drops rapidly, which is a result of the plastic deflection of the contact. In turn, the beam ends are loaded during the unloading process. For a loading force of F = 0.02kN, the beam loses contact with the rigid flat (Fig.10) in themiddle part, and thus only the ends of the beamare loaded. Finally, for a loading force of F =0.00kN the contact is completely unloaded, and the residual (plastic) deflection of the contact is observed.The reloadingprocess is very similar to unloading, but a higher force should be applied to cause the same contact deflections (Figs.8, 10, and 11). This model describes, for example, the contact pressure distribution wi- thin a screw joint. A similarmodel has been studied previously by researchers 520 R. Kostek Fig. 8. A force-displacement characteristic of the contact interface and the displacement y of the 11th node against the loading force F Fig. 9. The contact pressure distribution (a) and deflection of the contact (b) obtained for the initial loading Fig. 10. The contact pressure distribution (a) and deflection of the contact (b) obtained for the unloading The modelling of loading, unloading and reloading... 521 Fig. 11. The contact pressure distribution (a) and deflection of the contact (b) obtained for the reloading (Buczkowski and Kleiber, 2006; Marshall et al., 2006), but the issues of unlo- ading and reloading in the normal direction have not been considered. 3.3. Collision of two elastic bodies A collision between bodies is a classic issue in dynamics and is a problem that has been previously studied by many researchers. Interactions between bodies are usually describedwith simplemodels that reflect the restitution co- efficient, but a collision is a farmore sophisticated process. The kinetic energy of bodies can be dissipated by any of the following phenomena: plastic de- flection of contact, plastic deformation of bodies, hysteresis of the material (internal friction), propagation of cracks, and the friction force between the bodies.Moreover, the impact can excite vibrations of the bodies,which reduce the restitution coefficient (see e.g. Lifshitz and Kolsky, 1963; Oden and Mar- tins, 1985; Schiehlen et al., 2006;Wriggers, 2006). These phenomenamake the modelling of a collision difficult. In this Section, the collision of two elastic bodies made of steel is stu- died (Young’s modulus E = 2.1 · 1011Pa, Poisson’s ratio ν = 0.3, density ρ=7860kg/m3). The first body hits the second body with an initial velocity of Vx =0.5m/s, and the secondbody is fixedonone end (Fig.12).Bothbodies aremodelled with FEM constant strain triangular elements (CSTE), the con- stitutive relation of which is described byHooke’s law for plane stress. All the FEM elements are equilateral triangles, and their masses are lumped into the nodesof thefinite elementmesh(nodemasses) (Munjiza, 2004) (Fig.12).Thus, all acting forces are applied to the nodes. This approach leads to a number of second-order ordinary differential equations, which describe the dynamics of a point mass. This mass modelling technique provides an opportunity to 522 R. Kostek model the dynamics of two- and three-dimensional bodies (Stępniewski, 2007). The second body contains non-linear contact elements thatmodel the contact zone and contact forces (Fig.12). The normal forces are described with the contact model presented here, whereas the friction forces are described with the Coulomb model (coefficient of friction µ=0.1). Fig. 12. The system of two bodies (before being collided) modelled with FEM constant strain triangular elements and contact elements; me –mass of FEM element Each contact element contains twenty contact nodes. The contact nodes represent the contact of the roughness asperities, and so the contact forces are applied to the contact nodes (Fig.13c). Normal and friction forces are calculated using a special procedure. The first stage of the procedure verifies if the individual edges of the bodies are close to each other (Fig.13a). During this stage, the edges thatare close to eachother are selected, and theedges that are far from one another are eliminated. If the nodes of an edge and the nodes of the contact element are on the opposite sides of the line xmin, then they are not in contact (Fig.13a). The selection process is performed for four lines: xmin, xmax, ymin, ymax, which eliminates most of the edges. The elimination process is a type of the bounding box method (Munjiza, 2004) and is based on conditional expressions, which makes the procedure fast to execute. The edges are then translated and rotated (Munjiza, 2004). Then, the procedure calculates the contact deflection (penetration) for each contact node (Fig.13b). After this, the contact forces are calculated for each contact node (Fig.13c) The modelling of loading, unloading and reloading... 523 fromthevelocities andpositionsof thecontacting edges.Finally, the equivalent forces are applied to the nodes of the finite element mesh (Figs.11b,c). This description is abrief summaryof themain ideaof theprocedure.Themodelling of the contact by contact nodes can be used when a non-uniform mesh is applied or when two bodies have a contact element. Fig. 13. Stages of the procedure for calculating contact forces: selection of edges (a), calculation of contact deflection (b), and calculation of contact and nodal forces (c) The consecutive stages of the collision process are presented inFig.14.The initial velocity of the first body is Vx(t=0s) =0.5m/s and Vy(t=0s) =0.00m/s. A stage before the initial contact is shown in Fig.14a. Next, the magnitu- des of the Hubert-Mises-Hencky effective stress of the elements that are near the contact zone gradually increases because of the increasing contact for- ces (Figs.14a,b,c). The contact forces change the velocities of the first and second body at the same time (Figs.14a,b,c). A part of the initial kinetic energy of the first body is transferred to the second body, after which the first body rotates (Fig.14d). The contact zone has been loaded and unloaded during this period. The collision excites vibrations of the first body and the second body. The vibrations of the second body are presented in Figs.14d,e,f, and Fig.14g. The kinetic energy of the second body is converted to poten- tial energy, which causes the second body to deflect (Figs.14d,e,f). Thus, the velocity of the second body decreases, whereas the magnitude of the effec- tive stress increases. The potential energy is then transformed into kinetic energy (Figs.14f,g), which reflects the nature of vibrations. The two bodies are out of contact during this period. The average magnitudes of the velo- city components of the first body are V x(t=9.4·10−5s) = −1.25 ·10 −3m/s and V y(t=9.4·10−5s) =4.75·10 −2m/s. It shouldbementionedhere that at this stage, 524 R. Kostek Fig. 14. Consecutive stages of collision presented for the time t=0s (a) to t=3.850E-4s (l) The modelling of loading, unloading and reloading... 525 the y-component of the velocity is larger than the x-component. Finally, the secondbodyhits the first body (Figs. 12h and i), which changes the velocity of the first body (Figs.14h,i,j). The contact zone is reloaded andunloadedduring this secondcollision.Apartof thekinetic energyof the secondbody is transfer- red to the first body.After the second collision, the averagemagnitudes of the velocity components of the first body are V x(t=3.85·10−4s) =−3.30 ·10 −1m/s and V y(t=3.85·10−4s) =5.15 ·10 −2m/s. The vibrations of the bodies after colli- sions are presented in Figs.12k,l. Only a part of the kinetic energy is returned to the first body. The initial energy excited vibrations of the bodies and has beendissipated in the contact zone.Moreover, the y-componentof thevelocity appears as a result of friction phenomena. The results presented in this Section show the sophisticated nature of the collision process. Some simplifications used in the modelling of the collisions between the bodies, such as neglecting the deformations of the bodies or the application of simplified contactmodels, can lead to gross errors. The collision has been simulated and depicted using the proprietary programs FDEM RK and DEV KM (Kostek andMunjiza, 2009). The visualization of the presented collision is available on the Internet http://www.youtube.com/watch?v=ql28X42P5D0. 4. Conclusions A non-linear model of contact, which is dedicated especially to FEM, MBS and DEM, has been presented in this article. The model effectively demon- strates the contact characteristics, and it has a logical structure and a physical interpretation. Thephenomena that occur in the contact are non-linear, which makes the modelling of contact difficult. This model adequately describes the following characteristic phenomena: • the plastic deformation of a virgin contact, • the hysteresis loop of a preloaded contact, • the insensibility of the contact hysteresis to the frequency of loading. Themodel is effective and accurate and has been validated experimentally. Contact interfaces introduce plastic deflections, damping, hysteresis and memory of loading into many mechanical systems, which changes their stiff- ness, restitution coefficients and dynamical characteristics. Therefore, finding the relationshipbetween the loadinghistoryandthedeflectionhistory seems to 526 R. Kostek be important. Themodel presented in this paper can be applied to themodel- ling ofmechanical systems inwhich contact phenomenaplay a significant role, and can lead to the improved control of precision machines. Nowadays, a few microns of contact deflection are a significant factor for precision machining. Appendix This pseudo-code (Fig.15) presents an algorithmwhich calculates the contact deflection δ from the contact pressure p. In this case, the contact pressure is the input signal, whereas the contact deflection is the output signal. p, δ // input and output signals pmax, δep // variables which describe the inner state of the contact interface // for virgin contact they are zero δe = eep me; // elastic deflection if p>pmax then pmax = p else pmax = pmax δp = epp mp max; // plastic deflection if (δep/eep−e) 1/mep−e +(δep/eep−p) 1/mep−p >p> > (δep/eep−e) 1/mep−e − (δep/eep−p) 1/mep−p then dδep/dt=0; δep = δep; // elasto-plastic deflection, stick period if (δep/eep−e) 1/mep−e +(δep/eep−p) 1/mep−p

0; δep = eep−lp mep; // elasto-plastic deformation, loading curve if p< (δep/eep−e) 1/mep−e − (δep/eep−p) 1/mep−p then dδep/dt< 0; δep = eep−unlp mep; // elasto-plastic deformation, unloading curve δ= δp+δep+δe // contact deflection – output signal Fig. 15. The pseudo-code of the algorithmwhich calculates the contact deflection from the time history of contact pressure If the contact deflection δ is the input signal, then the contact pressure p is the output signal. An algorithm, which solves this issue, is more complica- ted than that above presented. This problem can be solved with a iteration procedure, e.g. false position method, or Newton-Raphson method. A special procedure dedicated to this issue can be applied as well. Acknowledgement I would like to take the opportunity to thank W. Woytek Tworzydlo, Ph.D., Director, Research and Engineering, Altair Engineering, Inc. Austin, TX, for the experimental results. The modelling of loading, unloading and reloading... 527 References 1. Archard J.F., 1957, Elastic deformation and the laws of friction,Proceedings of the Royal Society London, Ser. A, 243, 190-205 2. Buczkowski R., Kleiber M., 2006, Elasto-plastic statistical model of stron- gly anisotropic rough surfaces for finite element 3D-contact analysis,Computer Methods in Applied Mechanics and Engineering, 195, 37/40, 5141-5161 3. Cichowicz M., Nowicki B., 1984, Contact stiffness and experimental inve- stigations,XIII Sympozjum Tribologiczne, Częstochowa – Poraj, 316-321 4. Etsion I., Kligerman Y., Kadin Y., 2005, Unloading of an elastic-plastic loaded spherical contact, International Journal of Solids and Structures, 42, 13, 3716-3729 5. Goerke D., Willner K., 2008, Normal contact of fractal surfaces- Experimental and numerical investigations,Wear, 264, 7/8, 589-598 6. Greenwood J.A., Williamson J.B.P., 1966, The contact of nominally flat surfaces,Proceedings of the Royal Society London, Ser. A, 295, 300-319 7. Grudziński K., Jaroszewicz W., 2002,Posadowienie maszyn i urządzeń na podkładkach fundamentowych odlewanych z tworzywa EPY, ZAPOL 8. Grudziński K., Kostek R., 2007, An analysis of nonlinear normal con- tact microvibrations excited by a harmonic force, Nonlinear Dynamics, 50, 4, 809-815 9. GutowskiP., Skrodzewicz J., 2001,Experimental investigations of the pla- nar contact joints under dynamic loads, 5th Conference on Experimental Me- thods in Designing andMaintenanceMachines,WrocławUniversity of Techno- logy, Szklarska Poręba, 1, 341-352 10. Hess D.P., Soom A., 1991, Normal vibrations and friction under harmonic loads: Part II – Rough planar contacts,Transactions of the ASME, Journal of Tribology, 113, 1, 87-92 11. Hunt K.H., Crossley F.R.E., 1975, Coefficient of restitution interpreted as damping in vibroimpact, Transactions of the ASME, Journal of Applied Mechanics, 42, 2, 440-445 12. Jamari J., Schipper D.J., 2007, Plastic deformation and contact area of an elastic-plastic contact of ellipsoid bodies after unloading,Tribology Internatio- nal, 40, 8, 1311-1318 13. JonesR.E., 2004,Models for contact loadingandunloadingof a rough surface, International Journal of Engineering Science, 42, 17/18, 1931-1947 14. KadinY.,KligermanY., Etsion I., 2006,Unloading an elastic-plastic con- tact of rough surfaces, Journal of Mechanics and Physics of Solids, 54, 12, 2652-2674 528 R. Kostek 15. Kim J.Y., BaltazarA., Rokhlin S.I., 2004,Ultrasonic assessment of rough surface contact between solids from elastoplastic loading-unloading hysteresis cycle, Journal of the Mechanics and Physics of Solids, 52, 8, 1911-1934 16. KostekR., 2005, Investigationof thenormalcontactmicrovibrationsandtheir influence on the friction force reduction within the dynamical systems, Ph.D. Thesis, Technical University of Szczecin 17. Kostek R., Konowalski K., 2008,Amodel of contact between rough surfa- ces,Tribologia, 39, 5, 95-105 18. Kostek R., Munjiza A., 2009, Visualization of results received with the di- screte elementmethod,Computational Methods in Science andTechnology, 15, 2, 151-160 19. Kruggel-Emden H., Wirtz S., Scherer V., 2008, A study on tangen- tial force laws applicable to the discrete element method (DEM) for materials with visco-elastic or plastic behaviour, Chemical Engineering Science, 63, 6, 1523-1541 20. Li L.-Y., Gu J.-Z., 2009, An analytical solution for the unloading in sphe- rical indentation of elastic-plastic solids, International Journal of Engineering Science, 47, 3, 452-462 21. Lifshitz J.M., Kolsky K., 1963, Some experiments on anelastic rebound, Journal of the Mechanics and Physics of Solids, 12, 1, 35-43 22. Marshall M.B., Lewis R., Dwyer-Joyce R.S., 2006, Characterisation of contact pressure distribution in bolted joints, Strain, 42, 1, 31-43 23. Martins J., Oden J., Simoes F., 1990,A study of static and kinetic friction, International Journal of Engineering Science, 28, 1, 29-92 24. Michalczyk J., 2008, Phenomenon of force impulse restitution in collision modelling, Journal of Theoretical and Applied Mechanics, 46, 4, 897-908 25. MikicB.B., 1971,Analytical studies of contact of nominally flat surfaces effect of previous loading, Trans. ASME, Journal of Lubrification Technology, 20, 451-456 26. Mikic B.B., 1974, Thermal contact conductance; Theoretical considerations, International Journal of Heat and Mass Transfer, 17, 2, 205-214 27. MunjizaA., 2004,TheCombinedFinite-DiscreteElementMethod, JohnWiley & Sons Ltd. 28. NardinA., ZavariseG., SchreflerB.A., 2003,Modelling of cutting tools- soil interactions – Part I: Contact behaviour, Computational Mechanics, 31, 3/4, 327-339 29. Nicola L., Bower A.F., Kim K.-S., Needleman A., Van der Giessen E., 2007, Surface versus bulk nucleation of dislocations during contact, Journal of the Mechanics and Physics of Solids, 55, 6, 1120-1144 The modelling of loading, unloading and reloading... 529 30. Oden J.T., Martins J.A.C., 1985, Models and computational methods for dynamics friction phenomena, Computer Methods in Applied Mechanics and Engineering, 52, 527-634 31. PeiL.,HyunS.,MolinariJ.F.,RobbinsM.O., 2005,Finite elementmodel- ling of elasto-plastic contact between rough surfaces, Journal of the Mechanics and Physics of Solids, 53, 11, 2385-2409 32. Pombo J.C., Ambrósio J.A.C., 2008, Application of a wheel-rail contact model to railway dynamics in small radius curved tracks, Multibody System Dynamics, 19, 1/2, 91-114 33. Pullen J., Williamson J.B.P., 1972,On the plastic contact of rough surfa- ces,Proceedings of the Royal Society London, Ser. A, 327, 159-173 34. Schiehlen W., Seifried R., Eberhard P., 2006, Elastoplastic phenomena in multibody impact dynamics, Computer Methods in Applied Mechanics and Engineering, 195, 50/51, 6874-6890 35. Skrodzewicz J., 2003, Influence of the lubricating agent on the properties of contact joints, Journal of Theoretical and Applied Mechanics, 41, 1, 107-118 36. StępniewskiA., 2007,D’Alembert’s supplementedprincipleandNewton’sfive supplemented laws, International Journal of Pure and Applied Mathematics, 38, 3, 415-424 37. Stolarski T.A., 2000,Tribology in Machine Design, Elsevier 38. TworzydloW.W., CecotW., Oden J.T., YewC.H., 1998, Computatio- nalmicro- andmacroscopicmodels of contact and friction: formulation, appro- ach and applications,Wear, 220, 2, 113-140 39. Vu-Quoc L., Zhang X., 1999, An elastoplastic contact force-displacement model in the normal direction: displacement-driven version, Proceedings – Royal Society. Mathematical, Physical and Engineering Sciences, 455, 1991, 4013-4044 40. Woo K.L., Thomas T.R., 1980, Contact of rough surfaces: a review of expe- rimental work,Wear, 58, 2, 331-340 41. Wriggers P., 2006,Computational Contact Mechanics, 2nd Ed., Springer Modelowanie charakterystyki kontaktu uzyskanej podczas obciążania odciążania oraz powtórnego obciążania Streszczenie W artykule przedstawiono nieliniowy matematyczny model opisujący histerezę suchego kontaktu ciał chropowatych, który został obciążonyw kierunku normalnym. Kilka charakterystycznych zjawisk pojawia się podczas pierwszego obciążania, od- ciążania i powtórnego obciążania. Kontakt ciał jest najbardziej podatny podczas 530 R. Kostek pierwszego obciążania, a jego plastyczne odkształcenie jest znaczne. Natomiast cha- rakterystyka uzyskana podczas odciążania jest sztywniejsza, a odkształcenia mają głównie charakter sprężysty. Podobną charakterystykę uzyskuje się podczas ponow- nego obciążania, tak więc obserwowana jest pewna pętla histerezy. Zweryfikowany doświadczalnie model kontaktu, opisujący wyżej wymienione nieliniowe zjawiska, zo- stał przedstawionyw artykule.Model ten został opracowany zmyślą o:metodzie ele- mentów skończonych, metodzie elementów dyskretnych i układach wieloczłonowych. W praktyce obliczeniowej potrzebne są dokładne modele kontaktu, które nie wyma- gają jednak kosztownychobliczeń. Prezentację tegomodeluwzbogaconoprzykładami obliczeniowymi. Manuscript received August 11, 2011; accepted for print September 30, 2011