## 2 Materials and Methods
### Non-Hertzian Contact Model  
In order to obtain accurate pressure distributions, even in presence of a defect, a frictionless non-Hertzian numerical 3D surface contact model was used. This formulation was adopted since it is able to describe pressure distributions resulting from nonconstant curvatures and sudden geometry variations. Then, the discontinuities represented by defects are analyzed by means of this model. The contact conditions can be expressed in the so-called Hertz-Signorini-Moreau problem [61; 62; 63]. \[\mathbf{h}\geq\mathbf{0},\,\mathbf{p_{n}}\geq\mathbf{0},\,\mathbf{h}\cdot\mathbf{p_{n}}=\mathbf{0} \tag{5}\]  
The first condition enforces that no interpenetration can occur between the contact bodies, and therefore, the gap function \(\mathbf{h}\), which measures the distances between the surfaces, can only be positive or equal to \(\mathbf{0}\) in the contact area. In Equation (5), and in the following, bold lowercase variables represent vectors, while bold uppercase variables represent matrices. The second condition imposes that the contact is nonadhesive, and therefore, no tension force can be present in the contact area formulated from the normal stress:  
\[\sigma_{\mathbf{n}}=\mathbf{t}\cdot\mathbf{n} \tag{6}\]  
where \(\mathbf{t}\) is the traction force vector, \(\mathbf{n}\) is the normal direction to the surface and \(\mathbf{p_{n}}=-\sigma_{\mathbf{n}}\). The third condition enforces that the normal pressures can only be different from \(\mathbf{0}\) inside the contact area where \(\mathbf{h}=\mathbf{0}\) and null everywhere else. The gap function \(\mathbf{h}\) is expressed as  
\[\mathbf{h}=h_{0}+\mathbf{g}+\mathbf{\delta} \tag{7}\]  
where \(h_{0}\) is the indentation between the profiles imposed as a rigid body motion, \(\mathbf{g}\) is the initial separation of the contacting surfaces and represents its topography, while \(\mathbf{\delta}\) represents the elastic deformation of the surfaces due to the applied normal pressure \(\mathbf{p_{n}}\) and can be expressed as [64]  
\[\mathbf{\delta}=\mathbf{C}\cdot\mathbf{p_{n}} \tag{8}\]  
where \(\mathbf{C}\) is a matrix of the influence coefficients, which introduces the elasticity of the contacting surfaces. Its components \(C_{i,j}\) (\(i,j=0,1,\ldots,N\)) relate the displacement \(\delta_{i}\) at a point \(i\) due to the application of a unit pressure at point \(j\). The pressure distribution is estimated on the triangular grid-forming linearly varying pyramidal pressure elements, the definition of the influence coefficients \(C_{i,j}\) is detailed in [65; 66] and the solution procedure, which needs two iterative loops--one to eliminate tension forces and one to obtain the correct initial indentation \(h_{0}\)--is described in [54]. The computational model to which this contact algorithm will be applied is visible in Figure 2, alongside a detail of the representation of the geometry of the defect that will be introduced on the inner race for the reference size of 0.014 in (0.35 mm). Figure 2: (**a**) Computational model of the inner race and ball. (**b**) Detail of the damaged inner race for a defect size of 0.014 in (0.35 mm). Although only gravitational loads were supposed, since no direct information regarding the actual load applied to the damaged bearing was given and the details regarding the apparatus, the motor and the related weights, three load levels were hypothesized to have more data to extrapolate from. Those load levels were inferred from the bearing specifics and, indeed, from the manufacturer of the bearing; the fatigue load limit is given as \(P_{u}=335\) N, and therefore, the force \(F\) to be applied to the contact model was chosen as \(F_{1}=\frac{1}{S}P_{u}\), \(F_{2}=\frac{1}{4}P_{u}\) and \(F_{3}=\frac{1}{2}P_{u}\). The experimental activity was conducted in such a way that only the gravitational load acted on the bearings, but the information about weights and motor specifications were not available.