Multi-body Modelling of rolling eleMent bearings and perforMance evaluation with localised daMage

Rolling element bearings (REBs) are commonly used components in rotary machines. As their dynamic response has a great influence on the machine in which they are placed, they are key elements from the reliability point of view. To assure a satisfactory level of reliability during their useful life, an efficient maintenance approach such as condition-based maintenance (CBM) should be applied [16]. In this context, it is essential to detect, isolate and identify faults, as well as to determine the remaining useful life of the studied system; i.e. to carry out diagnosis and prognosis processes [2, 6, 16]. Physics-based models use physical laws and mathematical formulations to obtain the response of a system under certain operating conditions, making them an appropriate choice for such endeavours [14, 37]. These models can be used together with other methods, such as the data-driven approach or symbolic modelling, to create a hybrid model that aims to overcome the limitations of each method [9, 25]. In short, the development of physics-based models can facilitate the maintenance of REBs. Many researchers have considered the physics-based modelling of REBs. The majority of the literature on dynamic models of REBs focuses on ball bearings, with single-row deep-groove ball bearings receiving the most attention. Other dynamic models are related to other configurations of single-row ball bearings, such as angular-contact ball bearings [5, 15, 26, 42]. There are a few models for double-row ball bearings in the research literature [20, 46]. There are some models for roller bearings; these generally focus on cylindrical roller bearings [22, 28, 36], with a few some models proposed for tapered roller bearings [17]. Models for other configurations of roller bearings, such as double-row tapered roller bearings and double-row spherical roller bearings, are proposed by some authors as well [1, 3]. A review of the literature suggests there is a wider variety of physics-based models for ball bearings than roller bearings. Furthermore, litUrko LetUriondo oscar SALgAdo diego gALAr


Introduction
Rolling element bearings (REBs) are commonly used components in rotary machines.As their dynamic response has a great influence on the machine in which they are placed, they are key elements from the reliability point of view.To assure a satisfactory level of reliability during their useful life, an efficient maintenance approach such as condition-based maintenance (CBM) should be applied [16].In this context, it is essential to detect, isolate and identify faults, as well as to determine the remaining useful life of the studied system; i.e. to carry out diagnosis and prognosis processes [2,6,16].Physics-based models use physical laws and mathematical formulations to obtain the response of a system under certain operating conditions, making them an appropriate choice for such endeavours [14,37].These models can be used together with other methods, such as the data-driven approach or symbolic modelling, to create a hybrid model that aims to over-come the limitations of each method [9,25].In short, the development of physics-based models can facilitate the maintenance of REBs.
Many researchers have considered the physics-based modelling of REBs.The majority of the literature on dynamic models of REBs focuses on ball bearings, with single-row deep-groove ball bearings receiving the most attention.Other dynamic models are related to other configurations of single-row ball bearings, such as angular-contact ball bearings [5,15,26,42].There are a few models for double-row ball bearings in the research literature [20,46].There are some models for roller bearings; these generally focus on cylindrical roller bearings [22,28,36], with a few some models proposed for tapered roller bearings [17].Models for other configurations of roller bearings, such as double-row tapered roller bearings and double-row spherical roller bearings, are proposed by some authors as well [1,3].

sciENcE aNd tEchNology
tle research has been conducted on some REB configurations, including axial ball bearings, four-point contact ball bearings, double-row deepgroove ball bearings, slewing bearings, thrust roller bearings, etcetera.
In their study of single-row deep-groove ball bearings, Xiangyang and Wanqiang [43] used a two degrees-of-freedom (DOF) model considering the vertical and horizontal motion of the inner ring.Purohit and Purohit [30] took a similar approach, studying the effect of the number of balls and the preload on the frequency response of the REB.They concluded that an increase in the number of balls implies an increase of system stiffness and a reduction of the vibration amplitude.Kappaganthu and Nataraj [18] also studied motion in the plane of the REB, but they took into account the effect of internal clearance.
Focusing on the effect of defects, Qiu et al. [31] used a single DOF model to represent the dynamics of a REB in which both the system stiffness and damping are defined as the sum of undamaged and damaged parameters.The authors used different approaches to establish the evolution of damage, such as the linear damage rule, damage curve approach and the double linear damage approach.Rafsanjani et al. [32] reproduced the transient force that occurs when there is contact with a defective surface by means of a series of impulses repeated with the characteristic frequencies of the elements of the REB, whether a ring or a rolling element (RE).Patil et al. [29] and Kulkarni and Sahasrabudhe [21] introduced defects as geometric changes in the races using circumferential half sinusoidal waves and cubic Hermite splines, respectively.Pandya et al. [27] combined localised defects on the races and the REs and studied the dynamic response of a REB in those conditions.
The aforementioned models only take into account the motion in the plane of the REB.Other models, such as that used by Changqing and Qingyu [4], describe the three-dimensional (3D) motion using a five DOF model, defining three translation motions and two angular displacements.The mathematical modelling of the waviness is also considered to represent the geometrical imperfections.Zhang et al. [45] proposed a model for a rotor-bearing system using a five DOF approach and taking into account the effect of the lubrication by using an elastohydrodynamic lubrication (EHL) model in each ballrace contact.Sopanen and Mikkola [39] went beyond this approach and included both waviness and localised defects in their proposed model.Sawalhi and Randall [34] presented a bearing-pedestal model in which two DOF are used for the inner ring motion, two DOF for the pedestal motion and an extra DOF to simulate a high frequency resonant response.The authors included the effect of the slippage and localised faults in the races and the balls.Some authors have studied the motion of the shaft and the dynamics of the balls, defining one DOF [13,41] or two DOF [44] for each ball, but they have not considered the dynamics of the cage.In general, previous work is based on the assumption of the balls as equispaced.The motion of the cage is studied by a few authors; among these, Meeks and Tran [24] considered the dynamics of that component by using six DOF.
Dynamic models of single-row deep-groove ball bearings based on finite element modelling are available, such as those by Kiral and Karagüelle [19] and Liu et al. [23].In the latter work, the authors used a model consisting of 641,450 elements and 719,654 nodes with three DOF for each node, considering the cage to be rigid and the other components elastic.
As for other configurations of ball bearings, Wang et al. [42] used a three DOF model to determine the 3D motion of an angular contact ball bearing.They studied the radial, axial and angular displacements when a combination of radial, axial and moment loads are applied.Cui and Zheng [5] also considered this 3D motion but used a five DOF model to model a rotor-bearing system.Jain and Hunt [15] introduced a different approach; they calculated the friction forces between the balls and the races in each contact using the EHL theory.Although they did not consider the dynamics of the cage, they inserted springs of very high stiffness between the balls.Thus, the balls were not assumed to be equidistant.Niu et al. [26] introduced not only the lubrication tractive forces but also two translational DOF for the pedestal motion.Kogan et al. [20] combined the model of an angular contact ball bearing to obtain a model for a duplex angular contact ball bearing.Another model for a double-row ball bearing was offered by Zhuo et al. [46].The authors proposed a three DOF model for self-aligning bearings; they studied the effect of the applied loads, the internal clearance, the waviness and the number of balls on the dynamic response of the REB.
In work on roller bearings, Shao et al. [36] proposed a two DOF model for cylindrical roller bearings in which the time-varying deflection and stiffness are calculated when there is a localised defect in a race.Instead of studying the motion in the plane of the REB, the three DOF model of Patel and Upadhyay [28] considered 3D motion.These researchers also studied the effect of localised effects, but neglecting the edge stresses.Leblanc et al. [22] proposed a planar model for cylindrical roller bearings in which the cage is also modelled.They defined the normal forces as well as the friction torque in the contact between the rollers and the cage.This model includes the EHL theory for roller-race contacts.Kabus et al. [17] proposed a model for tapered roller bearings using the elastic half-space theory to define the contact pressure between the rollers and rings instead of the classical Hertzian contact theory.Models for other configurations of roller bearings, such as double-row tapered roller bearings and double-row spherical roller bearings, have been proposed by Bercea et al. [1] and Cao and Xiao [3], respectively.
The objective of this paper is to present a general approach to the physics-based modelling of REBs in different configurations, following a multi-body strategy.The model represents the characteristics of the parts that form the REBs as well as the way in which contacts between those parts are produced.It can obtain the dynamics of REBs of any kind of RE and in any configuration, based on the reusability of models.As different operating conditions affect both the shaft speed and applied loads, the model has the ability to reproduce either stationary or non-stationary operating conditions.The model is implemented using Modelica ® modelling language.
The paper is structured as follows: section 2 presents the basic concepts of the proposed physics-based model; section 3 shows how to create a model for any kind of REB based on the basics of the model; section 4 presents two case studies; section 5 provides a conclusion.

Definition of the basic concepts of the model
The multi-body model proposed in this paper is composed of three kinds of parts: the rings, REs and the cage.Each is assumed to behave as a rigid solid with point mass moving in the space, with the exception of the zones where there is contact between parts.Hence, each part is defined by six DOF: three for the translational motion and three for the rotational motion.These six DOF are defined by two vectors: r = [r x r y r z ] T and q = [q 0 q 1 q 2 q 3 ] T , where r is the linear position vector from the centre of a fixed reference system to the gravity centre (GC) of the body, and q represents the Euler parameters for the bodyfixed reference system [35].Only three of the four parameters are independent, as the condition of q being a unitary vector is imposed.
The dynamics of the aforementioned parts is affected by the contact loads and the loads caused by the restrictions.This section explains the physics behind these loads and the way damage is considered.

Contact loads between bodies
The contact between the different elements that form a REB causes forces and moments.Thus, each RE receives loads because of its contact with the rings and the cage, each ring receives loads because of its contact with each RE and the cage, and the cage receives loads from its contact with each RE and each ring.

Contact between the rolling elements and the rings
The contact force caused by the contact between the REs and the rings can be divided into two forces depending on the direction in which they are applied.One force is normal to the mating surface where the contact is produced; the other is on the contact plane.The directions of these forces are given by the unitary normal vector n and the unitary tangential vector t .
The normal force is calculated as the sum of elastic and dissipative components [7].The elastic force is determined using Hertz theory [12], whereby deformations occur in the elastic range and the dimensions of the contact area are small compared to the radii of curvature of the bodies under load.The dissipative effect is modelled as proposed by Flores et al. [7].The effect of the damping is introduced as a function of the relative normal penetration velocity δ  .Thus, the normal force f N is expressed as: (1 ) 4 where K n is the contact stiffness, δ the contact deformation, n the loaddeflection factor, which is equal to 3/2 for ball bearings and to 10/9 for any kind of roller bearings, ce the restitution coefficient, and δ  the initial normal impact velocity where the contact is produced.The expression to calculate K n is given as: where K i and K o are the stiffness of the contacts between the REs and the inner and outer rings, respectively.These two values are functions of the geometrical and material properties of the bodies, and the procedure to calculate them is explained by Harris and Kotzalas [12].The deformation is equal to zero if there is no contact and is calculated as the distance between the nearest surfaces of both the RE and the raceway of the ring, as expressed in the following equation: where r ψ is the distance between the GC of the RE and the raceway curvature centre in the radial plane, D w is the diameter of the RE, and r i,o is the raceway curvature radius of the ring.
Tangential force takes into account the effect of the lubricant and the rolling/sliding motion of the REs.According to the EHL theory, the shear stress τ in the lubricant film is defined as a function of the lubricant viscosity η and the strain rate γ  , i.e. the sliding speed between the two surfaces divided by the local film thickness.The integration of this shear stress in the contact area gives the friction force f T .In many cases, the friction force is obtained as the modulus of the normal force multiplied by a friction coefficient μ, i.e., as the shear stress τ, a function of η and γ  [40].Thus, In the proposed model, different expressions for the friction coefficient μ are used depending on the kind of RE.In the case of ball bearings it is obtained as [10,40] where τ e is the value of the shear stress at which shear thinning starts to become significant, p is the mean pressure in the contact, η 0 is the dynamic viscosity of the lubricant at atmospheric pressure, α is an average pressure-viscosity coefficient, U is the mean of the rolling speeds of the two surfaces with respect to the contact, SRR is the slide-roll-ratio [40], and h c is the central film thickness.Values for τ e and α depending on the mean film temperature and the mean contact pressure are given by Spikes [40].Hamrock and Dowson [11] offer an expression for the cental film thickness for elliptical contacts.
For roller bearings, another approach is followed when defining the value of μ [10], given as: where A, B, C and D are coefficients for a particular lubricant [10] and Δu is the sliding velocity of the surfaces.
The sum of the normal and tangential forces, leading to the total RE -ring contact force, creates a moment calculated as the cross product of the vector from the GC of each part to the contact points and the contact force.

Other contacts
The cage is considered as a body with Z holes for avoiding contact between the REs.These holes have a cylindrical form for ball bearings and a prism form for roller bearings.The contact loads between the REs and the cage holes are calculated using Hertz theory as explained above.
The contact between the rings and the cage is simplified by considering that a nonlinear spring connects the centre of the cage and the centre of the rings.This nonlinear spring is formed as a linear spring with a gap in which there is no contact force.Equations ( 7) and (8) show the relation between the radial distance between the centres r r and the radial component of the contact force f r , as well as the relation between the axial distance between the centres r a and the axial component of the contact force f a .The parameters defining those relations are the clearances c r and c a and the stiffness coefficients k r and k a for the radial and the axial force, respectively.As there is not considered to be any rotational stiffness in the contact between the rings and the cage, the moments are defined as equal to zero.

Housing restrictions and operating conditions
The housing is considered as either completely rigid or in a wide range of stiffness.The former results in the elimination of the six DOF of that ring; the latter implies the need to define a symmetric stiffness sciENcE aNd tEchNology matrix K h of order six with the aim of obtaining the loads because of the housing restrictions f h and m h , which are given as: where r r is the position vector of the ring and ϕ is the rotation vector of the ring, which is related to its Euler parameters [35].
The stiffness matrix K h is a function of the linear and angular position of the non-rotating ring and a function of its linear and angular velocity.The matrix can be calculated by estimating the coefficients or by reducing a finite element model of the housing to a matrix of order six.It should be noted that the rotating ring is considered to be rigidly fixed to the shaft.
The operating conditions are introduced as restrictions of the motion, in such a way that the rotary DOF is eliminated in order to impose a predefined speed on the rotary ring.It should be highlighted that either stationary or non-stationary operating conditions can be introduced.The applied loads are also given as a predefined profile function of time.They can be a combination of radial forces, axial forces and tilting moments and are added to the vectors f h and m h .

Introducing damage to the parts
Local damage is defined in the different parts of the REB as geometrical changes of their surface in one or more locations.The surface modifying profile function of the angular position directly affects the geometric properties of the contact between the elements.The deformation in a damaged case δ d is calculated as the subtraction or the sum of the deformation in the healthy case and the depth of the damage in the point at which the contact is produced, h, as expressed in the following equation: Equation ( 10) depends on whether there is a lack of material in a zone or some particle is fixed to any of the parts.Faulty situations such as spalling and wear debris in the contacts can be modelled using this approach.

Construction of the rolling element bearing model
Now that the general background has been established, this section explains the proposed approach.The first step is the creation of a new library of models, as shown in Table 1.These body and contact models are based on the use of attributes to simulate any kind of REB, as explained in Sections 3.1 and 3.2.
A model for a single-row REB is constructed by combining the appropriate models, as seen in Fig. 1.The schematic includes the bodies that compose a REB (two rings, Z REs and a cage), the contact models between the different elements, and the restriction models for the housing and the operating conditions.The lines between the boxes represent the connections between the body, contact and restriction models.The models of the bodies share their linear and angular position; these are used by the contact models to determine the existence of contact between the bodies and to establish the contact loads in an acausal way.The restriction models build conditions in the position of the bodies or in the loads supported by them according to Section 2.2.The motion of any part can be defined using the Newton-Euler equations.
With the aim of creating the model for REBs combining the different parts explained previously, this work uses Modelica ® modelling language.Modelica ® takes advantage of acausal modelling, its multi-domain modelling capability, its object-oriented approach and the ability to create and connect components of models to simulate dynamic systems represented by a set of differential algebraic equations (DAEs) [8].
As there are many different needs in industry with respect to rotary machines, many kinds of REBs have appeared over the years.This work seeks to model any kind of REB.To achieve this goal, the following issues are covered, based on the scheme of

Adapting the model for any kind of rolling element
One of the main features of a REB is the kind of RE used.REBs are classified as ball bearings or roller bearings; in addition, roller bearings can be composed of cylindrical, needle, tapered or spherical rollers.This paper defines two kinds of REs, as shown in Fig. 2. A sphere of diameter D w is defined for ball bearings (see Fig. 2a) and a general roller with mean diameter D w and length L w is defined for roller bearings (see Fig. 2b).
Despite their differences, the four kinds of rollers found in offthe-shelf REBs can be easily modelled using a roller model and vary-   2. The roller model only supports the combinations shown in Table 2; other combinations are not taken into account because it does not make sense.It should also be noted that the relationship between the length and the mean diameter is the only difference considered between the models of cylindrical rollers and needle rollers.Thus, the length of a needle roller is considered to be longer than four times its diameter.

Adapting the model for any internal configuration
Another important factor is the internal configuration of the REBs.Regarding this issue, many configurations can be found, that go from radial to thrust or axial configurations, through angular con-tact and four-point contact configurations in the case of ball bearings.
To obtain a model for a ring able to cover this wide range of configurations, a geometry of reference for rings is defined as shown in Fig. 3a.Rings for radial REBs, angular contact REBs and thrust REBs can be obtained by restricting the zone in which the contact between a RE and the ring is produced.Fig. 3b, 3c and 3d show, respectively, how to obtain outer rings for the aforementioned configurations based on the reference geometry.This restriction is given by the contact angle allowed to support each REB.Table 3 shows the ranges for the values of β 1 and β 2 , the initial and final possible contact angles, as shown in Fig. 3.The case of four-point contact REBs is constructed as the change in the number of curvatures to two, in such a way that the centres of curvature are displaced axially a distance L 4 , as shown in Fig. 4. Thus, the contact angles for each curvature are 0º < β 1 < β 2 < 90º for the raceway in the left and −90º < β 1 < β 2 < 0º for the raceway in the right.

Adapting the model for any number of rows
Once the kind of RE and the internal configuration of the REB are selected, the model for a single-row REB can be obtained.A multi-row REB is modelled by taking advantage of the reusability of the model.In general, a p-row REB is modelled taking p single-row REB models.To assure the p outer rings move in time with each other, they are connected by models containing stiffness matrices of order six, as shown in Fig. 5.The same procedure is used for the p inner rings and the p cages.These matrices relate the linear and angular position of the parts with the joint forces needed to assure the synchronous movement.The values of the matrices K o , K i and K g indicate whether the synchronisation between the parts is rigid or flexible.

Results and discussion
To show the capabilities of the proposed approach for the physics-based modelling of REBs, the paper employs two case studies.First, it gives results for the simulation of a single-row deep-groove ball bearing; more specifically, it simulates the REXNORD ER16K bearing.Second, it shows the response of a NJ 305 cylindrical roller bearing.The dimensions of the two REBs are shown in Table 4.
Both cases consider three different scenarios for the health condition of the REB: REB without damage, local damage on the outer ring and local damage on the inner ring.In the last two cases, damage is introduced as changes in the surfaces of the damaged bodies.As Fig. 6 shows, the damaged surfaces are defined by means of two smooth curves.These cubic sciENcE aNd tEchNology curves are defined by the angular length of the damage, ψ, and the damage depth, h.The properties of the material (steel) are the same for both case studies: modulus of elasticity E of 207 GPa, Poisson's ratio ν of 0.3 and density ρ of 7830 kg/m 3 .

First case study
The first case study considers a constant rotational speed of 20.9 rad/s (200 rpm) applied to the inner ring and a radial force of 1600 N. The selected values for ψ and h are 0.0847 rad and 0.13 mm, respectively, for both the damage in the outer ring and in the inner ring.
The model provides the time evolution of the different elements of the REB.Fig. 7 shows the vertical acceleration of the inner ring in the aforementioned three health conditions.As the figure indicates, the vibration of the healthy case is lower than that of the two damaged cases because of the presence of the damage in the surface of the rings.
To experimentally validate the results obtained by the model, the study uses a commercial test rig manufactured by SpectraQuest Inc., Gearbox Prognostics Simulator (see Fig. 8).An induction motor drives a motor in which a tachometer, a torquemeter and an encoder are mounted.In this case, one gearbox is monitored, more specifically, a REB in the intermediate shaft.There is another gearbox, as well as a motor that applies some load.The control guarantees that the speed of the drive motor and the load applied by the load motor are the desired ones.The vibration near the housing of the monitored REB is acquired using a triaxial PCB 356A17 accelerometer using a sample frequency of 10240 Hz.Fig. 9 shows the spectral content of the results obtained from the model and from the measurements taken in the test rig in both healthy and damaged conditions.For damaged conditions, localised damage in the outer ring is considered.In this case, periodic peaks are   found at a theoretical frequency called ballpass frequency of outer race (BPFO) [33], given as: where n i is the rotational speed of the shaft expressed in revolutions per second, D pw is the pitch diameter of the REB, and β is the contact angle.Fig. 9c shows the signal has some impulsiveness at the BPFO and at its next harmonics; the healthy case (Fig. 9a) shows a flatter spectrum with no peak at the aforementioned fault frequency.
A more complex signal can be seen when the fast Fourier transform is applied to the vibration measurements taken in the test rig, largely because of the different elements of the rig.In the results obtained from both the healthy and the damaged cases (see Fig. 9b and 9d) impulsiveness is found at the rotating frequency of the input shaft and the intermediate shaft of the monitoring gearbox (8.33 Hz and 3.33 Hz, respectively, shown with dashed lines) as well as some of their higher harmonics (6.66 Hz, 16.66 Hz, 25 Hz and 33.33 Hz, shown with dotted lines).The main difference between the healthy and the damaged cases is the presence of the BPFO (shown with a solid line) and higher harmonics (shown with dotted lines) in the latter because of the existence of damage in the outer ring.
The values of the peaks obtained in the four graphs shown in Fig. 9 are summarised in Table 5.As the table shows, the differences between the results are less than 0.85 %.This difference can be caused because of slip at the contact, as the theoretical value does not consider it.
Despite the lack of experimental results for the case with damage in the inner ring, the results obtained from the frequency-domain analysis are compared to the theoretical ones.Fig. 10 shows the spectral content of the acceleration of the inner ring with this kind of damage.The spectrum follows the expected pattern [38]; i.e., a peak is found at the theoretical frequency called ballpass frequency of inner race (BPFI) [33] expressed as: This is shown by a solid line.Impulsiveness is also found at the second harmonic of BPFI (36.13 Hz, shown with a dotted line), at the sidebands of BPFI and its second harmonic (11.23 Hz, 14.65 Hz, 21.48 Hz, 24.66 Hz, 29.54 Hz, 32.71 Hz, 39.31 Hz and 42.72 Hz, shown by dash-dot lines), at the rotating frequency (3.42 Hz, represented by a dashed line) and at the second harmonic of the rotating frequency (6.59 Hz, shown by a dotted line).In other words, the results obtained with the proposed model are near to the theoretical ones, with a difference up to 2.6 %.Data related to the BPFI and its second harmonic are summarised in Table 5.
The theoretical rotational speed of the cage ω c can be calculated as 2•π• FTF, where FTF is the fundamental train frequency [33], given as: The theoretical value of the rotational speed of the cage and the values obtained in the model for the different healthy scenarios are shown in Fig. 11 and summarised in Table 6 using the mean values of  the time signals.The differences are less than 0.35 % and are given by the loss of contact between the REs and the raceways when the REs reach the damaged zone; this produces a disturbance in the transmitted load and, therefore, in the friction forces producing the rotation of the balls involved in the motion of the cage.

Second case study
This case study considers the response of a cylindrical roller bearing when a radial constant force of 75 N and a constant shaft speed of 104.72 rad/s (1000 rpm) are applied to the inner ring.The selected values for ψ and h are 0.0093 rad and 0.5 mm, respectively, for the damage in the outer and inner ring, according to the tests done by Patel et al. [28].
Fig. 12 shows the results obtained from simulations of the vertical acceleration of the inner ring in the three health conditions.In order to experimentally validate the results of the model, results are compared to those obtained by Patel et al. [28] and to experimental data in the same work.Frequency-domain analysis is applied to the accelerations signals obtained from the proposed model for the two damaged cases; the results appear in Fig. 13.For the healthy case, an almost flat response is found; for the case with damage in the outer ring, impulsiveness is found at BPFO (solid line) and at the next two harmonics of this frequency (dotted lines); for the case with damage in the inner ring, peaks are found at the shaft frequency (dashed line), its higher harmonics (dotted lines), BPFI (solid line), two sidebands (dash-dot lines), and the second harmonic of BPFI (dotted line).The frequency values at which the characteristic fault frequencies and their corresponding second harmonics are found are summarised in Table 7; the table also contains the theoretical results and the simulation and experimental results of Patel et al. [28].The results obtained by the proposed model have differences of less than 1.4 % from the theoretical results, less than 2.3 % from the results of the model of Patel et al. [28] and less than 2.8 % from the experimental results shown in the same work.As Fig. 14 shows, the speed of the cage closely follows the theoretical values, with a mean difference of less than 1.2 %.The values of the cage speed are summarised in Table 8.Two case studies are analysed using the proposed approach.The response of a single-row deep-groove ball bearing and a cylindrical roller bearing is simulated.The vibratory response is indicative of the presence of damage in the REB and gives information on the kind of damage.The spectral content of the simulated results agree well with theoretical results, experimental data and results from other literature; hence, the model is validated.

Conclusions
The use of the proposed methodology for the physics-based modelling of REBs will help in future research into condition-based maintenance.Using the proposed model to carry out simulations covering a range of operating and health conditions will help to create a map of situations that may arise; a classification system can be set up, and the state of an operating REB, as well as its evolution, can be determined.

Fig. 1 :
How to model REBs using different REs.• How to model REBs in different internal configurations.• How to model REBs with different numbers of rows.•

Fig. 1 .
Fig. 1.Scheme of the model for a single-row REB

Fig. 7 .
Fig. 7. Vertical acceleration results obtained with the model in the first case study: (a) healthy case, (b) case with damage in the outer ring, (c) case with damage in the inner ring

Fig. 9 .
Fig. 9. Fast Fourier transform of the experimental and model results in the first case study: (a) model result for the healthy case, (b) experimental result for the healthy case, (c) model result for the case with damage in the outer ring, (d) experimental result for the case with damage in the outer ring This work proposes a methodology for the physics-based modelling of REBs.The objective is to represent the dynamics of different kinds of REBs and give detailed information about the response of a REB in different operating and health conditions using a multi-body model.To this end, the work takes advantage of the reusability of the different components of the model, generalised to represent as many cases as possible.Examples of this approach include modelling the REs, the different internal configurations and multi-row REBs.The model is able to consider the application of both stationary and non-stationary operating conditions to the rotary ring and the effect of damage considered as geometric changes in the surfaces of the bodies.It should be noted that other kinds

Fig. 12 .Fig. 11 .
Fig. 12. Vertical acceleration results obtained with the model in the second case study: (a) healthy case, (b) case with damage in the outer ring, (c) case with damage in the inner ring

Fig. 14 .
Fig. 14.Rotational speed of the cage in the different health conditions in the second case study

Table 1 .
Library of models

Table 2 .
Restrictions for the roller parameters in the different roller geometries

Table 3 .
Restrictions for the contact angles in the different internal configurations

Table 4 .
Dimensions of the analysed REBs inner raceway curvature radius, r i 4.1 mm outer raceway curvature radius, r o 4.1 mm -Fig.8. Test rig used for the experiments in the first case study sciENcE aNd tEchNology

Table 5 .
Comparison of BPFO and BPFI obtained from the model, using experiments and theoretical values in the first case study

Table 6 .
Comparison of cage speeds obtained from the model in the first case study in the different health scenarios and the theoretical

Table 7 .
Comparison between the fault frequencies obtained from the proposed model, from the work of Patel et al. [28] (both simulated and experimental) and the theoretical values in the second case study

Damage location Proposed model Theoretical Experimental data [28] Model of Patel et al. [28]
Hzof damage that do not imply local spalling are not included at this point.

Table 8 .
Comparison of cage speeds obtained from the model in the second case study in the different health scenarios and the theoretical value