# Modeling Impulsive Ball Mill Forces Effects on the Dynamic Behavior of a Single-Stage Gearbox

Gauthier Ngandu Kalala, Xavier Chiementin, Lanto Rasolofondraibe, Abir Boujelben, Bovic Kilundu

► To cite this version:

Gauthier Ngandu Kalala, Xavier Chiementin, Lanto Rasolofondraibe, Abir Boujelben, Bovic Kilundu. Modeling Impulsive Ball Mill Forces Effects on the Dynamic Behavior of a Single-Stage Gearbox. Machines, 2022, Machines Testing and Maintenance, 10 (4), pp.226. 10.3390/machines10040226. hal-03619674

HAL Id: hal-03619674

https://hal.univ-reims.fr/hal-03619674

Submitted on 25 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Modeling Impulsive Ball Mill Forces Effects on the Dynamic Behavior of a Single-Stage Gearbox

Gauthier Ngandu Kalala<fnref n="1" />,<fnref n="2" /> ID, Xavier Chiementin<fnref n="1" /> ID, Lanto Rasolofondraibe<fnref n="3" />,*, ID, Abir Boujelben<fnref n="1" /> and Bovic Kilundu<fnref n="4" />

**Abstract:** Gearboxes are frequently used in the mining industry, especially for power transmission between the electric drive and the ball mill; besides the extreme complexity of a ball mill gear transmission system, the fault diagnosis by vibration analysis can be easily distorted by the presence of impulsive noises due to the ball pulses on the mill shell. Although several works in the literature are related to the influence of an impulsive noise on the accuracy of the diagnosis, no dynamic model exists yet in the literature that can explain the influence of these forces on the dynamic behavior of gearboxes. This paper presents a new approach to determine the influence of the grinding forces in crack defects diagnosis. This approach is based on a hybrid numerical model of a 24-degree-of-freedom gearbox, simulating one gear train and two drive shafts. The impact forces of the mill drum are modelled by a discrete element method (DEM). The ball-filling rate ($F_r$), the mill speed ($N_r$), and the ball size ($D_b$) are considered to study this phenomenon. The simulations results show by a time series representation, fast Fourier transform, and short-time Fourier transform (STFT), that the acceleration is significantly affected by the presence of the grinding forces, developing an impulsive noise due to the impact of the balls governed by the studied parameters.

**Keywords:** ball mill; discrete element modeling; gear modeling; impulsive environment; vibration monitoring

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Ngandu Kalala, G.; Chiementin, X.; Rasolofondraibe, L.; Boujelben, A.; Kilundu, B. Modeling Impulsive Ball Mill Forces Effects on the Dynamic Behavior of a Single-Stage Gearbox. *Machines* **2022**, *10*, 226. https://doi.org/10.3390/machines10040226

Academic Editor: Davide Astolfi

Received: 1 March 2022

Accepted: 19 March 2022

Published: 23 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Ore treatment process is an essential step in the mining industry. The goal is to reduce it to a desired particle size by using a ball mill, to be then treated chemically or physically [1]. This process is most often associated with high operating costs. About 70% of these costs are attributed to the particle reduction range from 30–50 mm to 20–50 micron [2,3]. Research should be carried out to reduce these costs by exploring different solutions, for example, (i) optimizing the energy consumed by the ball mills [4,5], (ii) automation of the grinding process [6–8], or (iii) optimization of maintenance operations. The main purpose of industrial maintenance is to ensure the availability of the production tools. If well applied, it can also contribute efficiently to a reduction of the operating costs of the machines. Industrial maintenance has undergone a continuous evolution over the years, advancing from a corrective to a preventive approach. However, the drawback of this approach is the risk to intervene too late when a failure has already occurred before the next scheduled maintenance, or replacing parts too early, which might have been still functional before failure.

<footnote n="1">¹ Institute of Thermics, Mechanics and Material (ITheMM), University of Reims Champagne-Ardenne, Moulin de la Housse, 51687 Reims, France; gauthier.ngandu-kalala@univ-reims.fr (G.N.K.); xavier.chiementin@univ-reims.fr (X.C.); abir.boujelben@univ-reims.fr (A.B.)</footnote>

<footnote n="2">² Polytechnic Faculty, University of Lubumbashi, Kasapa, Lubumbashi 7110502, Congo</footnote>

<footnote n="3">³ Centre de Recherche en STIC (CReSTIC), University of Reims Champagne-Ardenne, Moulin de la Housse, 51687 Reims, France</footnote>

<footnote n="4">⁴ Mechanical Section, Higher Institute of Applied Technology (ISTA), Kinshasa 1021101, Congo; bovic.kulundu@ista.ac.cd</footnote>

<footnote>* Correspondence: lanto.rasolofondraibe@univ-reims.fr</footnote>
However, this could be corrected by a so-called conditional preventive approach according to which an intervention is decided based on the measured state of the equipment.

Vibration analysis has long been the main tool for monitoring rotating machines and is therefore highly represented in the industry. It operates based on the properties of impulsivity and periodicity of disturbances related to a localized fault. However, when applied to nonstationary operating conditions and impulsive environments such as in the mining industry, it shows a lack of precision to detect mechanical defects. The monitoring methods using vibration analysis developed in the context of rotating machines should be accurate in order to ensure a reliable diagnosis. To obtain satisfactory results, real conditions in which the monitored equipment operates must be considered. Through modeling and simulation, it is becoming increasingly possible to understand certain complex phenomena normally rather difficult to obtain directly in an industrial environment. This paper focuses on the vibration behavior of gearboxes driving ball mills in the mining industry, which often operate in a harsh environment due to shocks and collisions on the mill drum.

Research has already focused on these issues, particularly in the context of gear modeling. The first trend focuses on the consideration of varying operating conditions. This research seeks to understand, with the help of modeling, the consequences that fluctuations in loading and/or speed can have on the diagnosis of gear defects by vibration analysis. F. Chaari et al. [9] proposed a model for the variation of the gearing stiffness based on the mechanical characteristics of the drive motor and the loading conditions. Simulation results on a single-stage spur gear transmission led to frequency modulations in perfect agreement with experimental results. W. Bartemus and R. Zimroz [10] demonstrated through a concentrated mass modeling and an experimental investigation the sensitivity of a gear transmission in terms of vibration level and increasing loading conditions. They showed that a defective transmission is more sensitive to increased load by generating higher vibrations. N. Baydar and A. Ball [11] concluded that a spectral analysis cannot track the degradation of a gear's condition under fluctuating load conditions. They used the instantaneous power spectrum to detect local defects in the teeth of a gear under different load conditions. Farhat et al. [12] developed a numerical model of single-stage gears subjected to variable speed and load conditions and studied their influences in the presence of combined gear and bearing faults.

Other researchers have been interested in the environment, in which most rotating machines operate, mainly in the mining industry. The presence of shocks, impulses, and collisions of materials that characterize most mining processes is often highlighted. It has been observed that the vibration signatures collected in these processes are most often contaminated by the presence of noncyclic impulsive shocks at high energy. This phenomenon makes it difficult to extract useful information about the monitored machine faults. Stochastic signal modeling approaches are most often used for these types of processes. G. Yu and N. Shi [13], for example, proposed a new method of statistical modeling of gear defects based on an α-stable distribution. A. Wylomanska et al. [14] also proposed a method for separating the useful signal related to the bearing fault from a noncyclic and non-Gaussian impulsive noise, in the context of monitoring a raw material crusher. The proposed technique was applied to a simulated signal as well as a real signal. S. Schmidt et al. [15] suggested the synchronous median of the square of the envelope method, instead of the synchronous mean of the square of the envelope of a vibration signal, for fault diagnosis of gearboxes operating under variable conditions and in an impulsive environment. This study showed its effectiveness on simulated and real signals. Jacek w et al. [16] proposed to use nonnegative matrix factorization of spectrogram for separation cyclic and noncyclic impulsive components from a hammer crusher. Jacek w. [17] proposed the time-varying spectral kurtosis (TVSK) as a means of spectral kurtosis adaptation for frequency-band selection of bearing damage signals under variable operating conditions. In a corresponding trend, Hebda-Sobkowicz et al. [18] presented a new approach to local bearing damage detection in the presence of non-Gaussian impulsive noise based on conditional variance statistics to identify cyclic and noncyclic impulses. This method
allows to detect and extract cyclic-impulsive signal (damage in bearing) in the presence of
high amplitude noncyclic impulsive signal. Several other researchers have addressed the
issue in this sense [19–21].

Gearboxes driving ball mills in the mining industry are subject to impulsive forces created by the falling of the balls and ore against the mill drum. The vibratory characteristics of the gears leading to a fault detection can be affected by the presence of these forces/shocks, which makes it difficult to detect localized faults in the gearboxes. However, in the literature, no numerical model allows us so far to understand the vibratory behavior of gearboxes in an impulsive environment such as that generated by the ball mill. The lack of knowledge of the nature of the forces generated by the fall of the balls against the mill drum does not permit us to develop an efficient mathematical model that characterizes the distribution of these forces. For modeling gears in an impulsive environment, the common practice is to synthesize the signal, considering it as a combination of several components: a deterministic component which often characterizes the manifestation of the fault, a gaussian random component which models noise, and a non-Gaussian random component which models impulsive noise. The weakness of this method is that the modeled signal can only describe one type of defect at a time. In this paper, a new hybrid numerical approach is proposed, combining two types of models: (i) a classical model of a gear transmission (lumped parameter model and finite elements models for gear and shaft), (ii) a discrete element model of a ball mill. The second model simulates the internal behavior of the ball mill and provides the grinding forces that will be considered as external efforts in the model of the transmission. The objective is to show, using this model, the influence of impulsive grinding forces on the transmission vibratory behavior in the presence of defects. In addition to the introduction and the conclusion, the paper is organized into three main sections. Section 2 presents the general methodology of the hybrid approach giving details on the software and the design of the simulation to be performed. Section 3 focuses on the description of the gearbox model and the ball mill model separately. Section 4 presents results; first, the load profile and the shape of the grinding forces are presented, then the influences of its forces are studied in healthy and under crack fault conditions, respectively.

**2. Methods**

2.1. *Introduction*

In industrial reality, a ball mill is generally made of a cylindrical drum driven by an
electric motor via a gear reducer. In addition to the gear reducer, a ring gear is directly
connected to the mill, whose transmission is provided by an external pinion, thus yielding
an additional reduction stage. The model studied in this paper is a simplified version
in accordance with a laboratory setup in which the mill is mounted at the end of the
gearbox output shaft (Figure 1). During its rotation, the mill drum lifts the charge along
one side until it reaches the maximum drop height. From this position, the particles start
to fall independently of the drum movement, described by either a cascade or a cataract
movement depending on the rotation speed of the mill [3].

By falling against the drum, particles' impact forces are created and transmitted to the gearbox through the drive shaft. Despite the presence of a lining to absorb the impacts, the impact forces produced by the movement of the load inside the mill are present and their amplitude can increase with the wear of the lining. It is therefore necessary to study their influence on the dynamic behavior of the gearbox.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Ball mills driving system.</figcaption></figure>

## 2.2. Methodology

The methodology implemented in this work is presented in the (Figure 2). It is based on a modeling approach called hybrid, which consists of bringing together different types of models for each technological element involved in the system. Three technological elements were involved in this model. The drive shafts were modeled by a finite element method, using Timoshenko beam elements. Gears were described by a lumped parameter model and were located on their corresponding nodes on the two input and output shafts; the ball mill was represented by a discrete element method. The geometry of their shell, produced by computer-aided design software such as SolidWorks and/or FreeCAD, can be loaded from many commonly used formats such as a mesh produced by Gmsh software [22]. The discretized model was loaded using a Python program; spherical particles were filled in a ball mill and run using discrete elements modeling in Yade software. By solving Newton's second law, the trajectory of every particle was calculated, so the properties and behavior of the material under the influence of different external forces may be observed. By observing the trajectory and postprocessing of the collected data either in Yade or an external program such as ParaView, a deeper understanding of the microscopic processes may be obtained [23]. After obtaining impulsive forces and assembling each element, the calculated forces were applied to the 24-degree-of-freedom gearbox model and were solved using the Newmark time integration scheme in (MATLAB), to study the dynamic response.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Methodology.</figcaption></figure>
## 2.3. Simulation Design

During the operation of the mill, the rotational speed, filling rate, and ball size are operating parameters that can change during the process and influence in different ways the dynamics of the load inside the mill, and consequently, of the entire structure [4]. The mill speed is defined as the percentage of the critical speed $C_s$. The critical speed $C_s$ of the mill is the rotation speed at which the balls are subject to centrifugal force while adhering to the inner surface of the drum [4,24,25]. It can be calculated by the following Equation (1).

$$C_s = \frac{42.3}{\sqrt{D_M}} \qquad (1)$$

where $D_M$ is the average diameter of the mill. The critical speed of the mill simulated in this paper was 5.04 (rad/s). Four speed levels from 30% to 95% of this speed were considered. The filling rate $F_r$ is the fraction of the internal volume of the mill occupied by the circulating load. The four simulated filling rates ranged from 5% to 35% and the simulated mill was successively fed with balls of diameter ranging from 30 to 120 mm. The choice of these three parameters is justified by the fact that their influence has been noticed through industrial observations during the grinding process. For example, the reduction in ball size due to wear leads to an outward ejection of the balls. This phenomenon has a direct influence on the filling rate. It is often decided to increase the number of balls to achieve the filling rate required for grinding ores. Moreover, the rotation speed of the mill affects the type of ore grain size and varying this speed can often adjust the grain size at the outlet of the mill by. The grinding can therefore be carried out in cascade or cataract modes, depending on whether the rotation speed is close to or far below the critical speed. Table 1 shows how the simulations were designed according to the three operating parameters.

<table><caption>**Table 1.** Simulation design.</caption><thead><tr><th>Operating Conditions</th><th>Rotating Speed</th><th>Ball Filling Rate</th><th>Ball Sizes</th></tr></thead><tbody><tr><td>Fixed parameters:</td><td>F<sub>r</sub> = 30%, D<sub>b</sub> = 50 mm</td><td>N<sub>r</sub> = 0.6 C<sub>s</sub>, D<sub>b</sub> = 50 mm</td><td>F<sub>r</sub> = 30%, N<sub>r</sub> = 0.6 C<sub>s</sub></td></tr><tr><td>No.</td><td>% of Critical speed C<sub>s</sub></td><td>Filling rate, F<sub>r</sub></td><td>Diameter (mm), D<sub>b</sub></td></tr><tr><td>1</td><td>30%</td><td>5%</td><td>30</td></tr><tr><td>2</td><td>50%</td><td>15%</td><td>60</td></tr><tr><td>3</td><td>70%</td><td>25%</td><td>90</td></tr><tr><td>4</td><td>90%</td><td>35%</td><td>120</td></tr></tbody></table>

# 3. System Modeling

## 3.1. Modeling of the Gearbox

* The transmission system studied in this work was a single-stage gearbox, according to the bench SURVIB of the laboratory of mechanical engineering of the URCA. It consists of two transmission shafts, two gears, four bearings, and a housing. Some assumptions were made in the modeling of the gears:

* The influence of the housing was neglected.

* The shafts were assumed to have constant cross sections.

* The wheel and the pinion were rigid, the only deformation being at the contact point, and the connection was flexible, with a variable mesh stiffness.

* The contact between wheel and pinion was assumed to be frictionless.

### 3.1.1. Gear Element

The chosen gear model considered the work previously carried out by [26] in the context of wind turbine gearboxes. In his work, the developed model allowed us to connect the node $N_1$ of the wheel to the node $N_2$ of the pinion as described in (Figure 3b), by linking
through a stiffness matrix $K_{gear}$ the restoring forces transmitted to the displacements $X_{gear}$, (Figure 3a), as expressed by Equations (2) and (3).

$$ \{F_{gear}\} = [K_{gear}] \{X_{gear}\} \qquad (2) $$

and

$$ [K_{gear}] = k_m [G] \qquad (3) $$

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3. Dynamic model of the gearbox. (a) Description of gear connecting forces; (b) representation of gear nodes disposition.</figcaption></figure>

The stiffness matrix $[K_{eng}]$ of the gear element can be expressed as a product of a matrix $[G]$ including the elements describing the geometry of the gear (radii of the base circles, pressure angle $\alpha$, and an average mesh stiffness $k_m$, which can be estimated according to ISO 6336 [27]). The matrix $[G]$ is calculated by considering that the strain energy associated with this average stiffness is equal to:

$$ E_p = \frac{1}{2} k_m \cdot \Sigma^2 = \frac{1}{2} k_m (\{\sigma\}^T \cdot \{X_{gear}\})^2 = \frac{1}{2} k_m (\{X_{gear}\} [G] \{X_{gear}\}^T) \qquad (4) $$

$$ [G] = \sigma^T \sigma \qquad (5) $$

where $\sigma$ is the geometric vector of the gear transmission and defined by Equation (6):

$$ \sigma = \{ c \ 0 \ c \ -\lambda_1 \ \mu_1 s + d_1 c \ \lambda_1 c \ -c \ 0 \ -s \ \lambda_2 s \ -d_2 c \ -\lambda_2 c \} \qquad (6) $$

$\{X_{gear}\}$ is the generalized gear displacement vector defined by:

$$ \{X_{gear}\} = [x_p \ y_p \ \theta_p \ x_g \ y_g \ \theta_g]_{gear} \qquad (7) $$

The parameters involved in Equation (6) are described in Table 2 for a two-mesh configuration depending on the type of loading.

$R_p$ and $R_g$ are, respectively, the radii of the pinion and wheelbase circle; $\alpha$ and $\beta$, are, respectively, the pressure angle and the helix angle of the tooth in contact.

The modeling led to the determination of the mesh rigidity matrix as described in Appendix A. Compared to the work of [26] in which the meshing stiffness was assigned a constant, a variable meshing stiffness was considered in this work, in order to describe the manifestations of meshing defects. However, nonlinear bearing forces were not dealt with in this case. Due to the variation of the number of teeth in contact during meshing,
the duration of the contact depends on the contact ratio $\epsilon$ and the variation of the meshing stiffness can be modeled by a square function (Figure 4) K<sub>gear</sub>(t), according to the work done by [9,12]. It fluctuates around a mean value k<sub>m</sub> and can be described by Equation (8)

$$ K_{\text{gear}}(t) = \begin{cases} K_2 & 0 \le t \le (\epsilon - 1)T_m \\ K_1 & (\epsilon - 1)T_m \le t \le T_m \end{cases} \qquad (8) $$

where T<sub>m</sub> is the meshing period, ε is a contact ratio, and K<sub>1</sub> and K<sub>2</sub> are the minimum and maximum values corresponding to the contact stiffness when one or two teeth are in meshing. A crack defect can be introduced to the model, by cyclically weakening the mesh stiffness by a degradation rate calculated by Equation (9):

$$ D(t) = \alpha\delta \left[ 1 + \sum_n \frac{\sin(n\alpha\pi)}{n\alpha\pi} \cos(n\Omega_i t) \right] \qquad (9) $$

<table><caption>**Table 2.** Gear configurations.</caption><thead><tr><th>Parameters</th><th>Configuration 1</th><th>Configuration 2</th></tr></thead><tbody><tr><td>μ<sub>1</sub></td><td>-R<sub>p</sub> tan α</td><td>-R<sub>p</sub> tan α</td></tr><tr><td>μ<sub>2</sub></td><td>+R<sub>g</sub> tan α</td><td>+R<sub>g</sub> tan α</td></tr><tr><td>λ<sub>1</sub></td><td>R<sub>p</sub></td><td>-R<sub>p</sub></td></tr><tr><td>λ<sub>2</sub></td><td>-R<sub>g</sub></td><td>-R<sub>g</sub></td></tr><tr><td>c</td><td>cos β</td><td>cos β</td></tr><tr><td>s</td><td>sin β</td><td>sin β</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Fluctuation of meshing stiffness K<sub>gear</sub>(t).</figcaption></figure>

In the presence of this defect, the meshing rigidity is thus calculated by:

$$ K(t) = K(t)[1 - D(t)] \qquad (10) $$

with 0 < D(t) < 1, and Ω<sub>i</sub> = 2πf<sub>i</sub>, f<sub>i</sub>, representing the rotation frequency of the wheel affected by a crack. δ and α = 1/Z<sub>i</sub> are the amplitude of the function D(t), and the cyclic ratio, respectively. Table 3 gives the different parameters used for the gear simulation.

<table><caption>**Table 3.** Simulation parameters of the gearbox.</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Wheel</th><th>Pinion</th><th>Motor</th><th>Units</th></tr></thead><tbody><tr><td>Mass</td><td>M</td><td>2.5</td><td>1.8</td><td>8</td><td>Kg</td></tr><tr><td>Moment of inertia</td><td>J</td><td>0.0014</td><td>0.007</td><td>0.00125</td><td>k·gm<sup>2</sup></td></tr><tr><td>Module</td><td>m</td><td>2</td><td></td><td></td><td>mm</td></tr><tr><td>Angle of pression</td><td>α</td><td>20</td><td></td><td></td><td></td></tr><tr><td>Number of teeth</td><td>Z</td><td>59</td><td>66</td><td></td><td></td></tr><tr><td>Meshing stiffness</td><td>K<sub>m</sub></td><td>1.5 × 10<sup>8</sup></td><td></td><td></td><td></td></tr><tr><td>Base circle radius</td><td>R</td><td>59</td><td>66</td><td></td><td>mm</td></tr></tbody></table>
### 3.1.2. Shaft Element

The shaft of the gearbox was assumed to be flexible in bending and torsion and made of a homogeneous material, to have a constant cross section, and to have small displacement in bending. An energetic approach based on the Lagrange formalism was applied. The two shafts were discretized into 6 elements including 8 nodes, as shown in (Figure 5) to approximate the displacement field. Unlike the work carried out by [28], the 3-degree-of-freedom in Equation (11), corresponding to one rotation along the z-axis and two translations along x and y were retained in each node to describe the system. The nodal displacement vector {$\chi$} is defined in Equation (12):

$$ \{\chi\}_{\text{node}} = \{x, y, \theta\}_{\text{node}} \qquad (11) $$

<figure><img src="image_9.png" /><figcaption>Figure 5. Discretization of the shaft elements.</figcaption></figure>

The vector containing the displacement fields over an element is given by:

$$ \{\chi^e\} = \{x_1, y_1, \theta_{z1}, x_2, y_2, \theta_{z2}\} \qquad (12) $$

Thus, the model has a total of 24 degrees of freedom (Figure 5). The adequate choice of shape functions and application of Lagrange's equations led to the determination of elementary mass **Me** and rigidity **Ke** matrixes of the shaft element, as described in Appendix B. The motor and the load were accounted for by their respective masses and polar inertias. They were concentrated at node 1 of the input shaft and node 8 of the output shaft, respectively. The bearings, modeled by constant rigidities mounted in parallel, supported the shafts. The geometric and material properties of the shafts are reported in Table 4.

<table><caption>Table 4. Simulation parameters of shaft.</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Input Shaft</th><th>Output Shaft</th><th>Units</th></tr></thead><tbody><tr><td>Shaft length</td><td>L.</td><td>0.5</td><td>0.5</td><td>m</td></tr><tr><td>Wheel position</td><td>Lg</td><td>0.3</td><td>0.3</td><td>m</td></tr><tr><td>Shaft outer radius</td><td>r<sub>e</sub></td><td>0.018</td><td>0.018</td><td>m</td></tr><tr><td>Mass per unit volume</td><td>ρ</td><td>7800</td><td>7800</td><td>k·gm<sup>-3</sup></td></tr><tr><td>Young modulus</td><td>E</td><td>2.10<sup>11</sup></td><td>2.10<sup>11</sup></td><td>N·m<sup>-2</sup></td></tr></tbody></table>
### 3.1.3. Discrete Element Modeling of the Ball Mill

#### 3.1.3.1. Principle

The discrete element method is used to simulate the motion and interactions of discrete particles in a dynamic environment such as a ball mill. This method uses contact models for an accurate characterization of the contact forces between the particles during their motion. In this paper, a process based on the discrete element method was developed and implemented in Yade in order to determine the motion of the charge and to calculate the interaction forces between the ball and the mill drum. The numerical model of the mill had therefore two main components: first, a geometric component that represents the geometry of the surface with which the balls come into contact (Figure 6a); second, a contact model that describes the interactions between the modeled particles and a drum as well between the different particles (Figure 6b). Three different elements were considered in the simulation: spherical elements representing the balls, triangulated surfaces representing shell liner and planes representing end walls.

<figure><img src="image_10.png" /><figcaption>Figure 6. Discrete element modeling of ball mill: (a) load movement profile inside de ball mill; (b) contact model.</figcaption></figure>

Considering a set of particles whose motion is described by Newton's second law, the linear contact model shown in Figure 6b models the forces acting on each particle. In ball mills modeling, these particles are usually spheres representing the balls. The following properties characterized each sphere j: (i) properties of states such as the position of the center $x_{bj}$, the linear speed $\dot{x}_{bj}$, the angular speed $\omega_{bj}$, the linear acceleration $\ddot{x}_{bj}$, the angular acceleration $\dot{\omega}_{bj}$, the forces applied to their centers $F_{bj}$, and the moment of force, i.e., the torque about the center $M_{bj}$; (ii) shape properties such as diameter $d_{bj}$ or the radius $r_{bj}$ of each particle; (iii) material properties such as density $\rho_{bj}$ and contact parameters. These spheres move in space according to Newton's second law (in translation and rotation) as follows:

$$ \left\{ \begin{aligned} m_{bj}\ddot{x}_{bj} &= F_{bj} \\ I_{bj}\dot{\omega}_{bj} &= M_{bj} \end{aligned} \right. \qquad (13) $$

where $m_{bj} = \frac{4\pi\rho_{bj}r_{bj}^3}{3}$, $I_{bi} = \frac{2m_{bj}r_{bj}^2}{5}$, are, respectively, the mass and the moment of inertia, and $F_{bj}$ represents the mechanical contact forces. For elastic particles, $F_{bj}$ was determined using the classical contact theory.
### 3.1.3.2. Model of Contact between Particles

The interactions between the balls and the mill wall can be dealt with by means of a suitable contact law. Both gravity and friction effects were considered. A linear contact law was used to determine the normal contact forces $F_n$, and tangential forces $F_s$ [29];

The normal force $F_n$ was calculated based on the relative displacement $\delta_n$, created during the interaction between the sphere and the contact surface. It was expressed as follows:

$$F_n = K_n \delta_n + \eta_n \frac{d\delta_n}{dt} \qquad (14)$$

where $K_n$ and $\eta_n$ are, respectively, the stiffness and the viscosity coefficients in the normal direction. The tangential force $F_s$ was defined by the following relation:

$$F_s = \min(\mu F_n, K_s \int v_s dt + \eta_s v_s) \qquad (15)$$

It consists of the frictional force $\mu F_n$, where $\mu$ is the friction coefficient acting between the particles, and the force due to the spring and the damping in the shear direction exceeding $\mu F_n$. $\eta_s$ represents the viscosity coefficient in the tangential direction and $v_s$ is the velocity in the same direction. Table 5 gives different parameters used in the ball mill simulation.

<table><caption>**Table 5.** Discrete element method simulation parameters.</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Radii of the balls</td><td>R<sub>i</sub></td><td>50</td><td>mm</td></tr><tr><td>Ball density</td><td>ρ<sub>i</sub></td><td>7640</td><td>Kg/m<sup>3</sup></td></tr><tr><td>Ball mill liner density</td><td>ρ<sub>b</sub></td><td>7640</td><td>Kg/m<sup>3</sup></td></tr><tr><td>Friction coefficient (balls–liners)</td><td>μ<sub>1</sub></td><td>0.30</td><td></td></tr><tr><td>Friction coefficient (ball–ball)</td><td>μ<sub>2</sub></td><td>0.30</td><td></td></tr><tr><td>Restitution coefficient</td><td>e<sub>1</sub></td><td>0.75</td><td></td></tr><tr><td>Restitution coefficient (balls–liners)</td><td>e<sub>2</sub></td><td>0.75</td><td></td></tr><tr><td>Mill diameter</td><td>D</td><td>0.80</td><td>m</td></tr><tr><td>Mill width</td><td>L</td><td>0.40</td><td>m</td></tr><tr><td>Normal rigidity</td><td>K<sub>n</sub></td><td>10<sup>6</sup></td><td>N/m</td></tr><tr><td>Tangential rigidity</td><td>K<sub>s</sub></td><td>10<sup>6</sup></td><td>N/m</td></tr><tr><td>Damping normale coefficient</td><td>η<sub>n</sub></td><td>10<sup>6</sup></td><td>Ns/m</td></tr><tr><td>Damping shear coefficient</td><td>η<sub>s</sub></td><td>10<sup>6</sup></td><td>Ns/m</td></tr></tbody></table>

The equations of motion described by Equation (13) were solved numerically, using an explicit integration scheme called the “Leapfrog” algorithm implemented in Yade software, and described by the following equation.

$$ \begin{cases} \dot{q}\left(t + \frac{\delta t}{2}\right) = \dot{q}\left(t - \frac{\delta t}{2}\right) + \ddot{q}(t)\delta t \\ q(t + \delta t) = q(t) + \dot{q}\left(t + \frac{\delta t}{2}\right) \end{cases} \qquad (16) $$

where $q$, $\delta t$, and $t$ are, respectively, the generalized position vector of the sphere mass centers, the constant time step, and the current time. This algorithm is conditionally stable. To ensure the stability of the algorithm, the time step should be less than the critical time step $\delta t_{cr}$ given by:

$$ \delta t_{cr} = \frac{2}{\omega_{max}} \qquad (17) $$

$\omega_{max}$ being the highest frequency in the system [30].
### 3.1.4. Equations of Motion

The above-described elements were assembled to obtain the finite element model. The latter was governed by a 24-equation system corresponding to the 24 degrees of freedom described in Section 3.1.2. Equation (18) represents the described model as a matrix equation.

$$[\mathbf{M}]\{\ddot{\mathbf{X}}\} + [\mathbf{C}]\{\dot{\mathbf{X}}\} + [\mathbf{K}]\{\mathbf{X}\} = \{\mathbf{F}_{\text{BM}}(\mathbf{t})\} + \{\mathbf{F}_{\text{ext}}(\mathbf{t})\} \quad (18)$$

where $\mathbf{F}_{\text{ext}}(\mathbf{t})$ is the vector of external forces applied to the gearbox, consisting of the motor torque $\mathbf{C}_{\text{M}}(\mathbf{t})$, and the resistive torque of the mill $\mathbf{C}_{\text{R}}(\mathbf{t})$. $\mathbf{F}_{\text{BM}}(\mathbf{t})$ represents the vector of forces exerted by the mill.

$$\mathbf{F}_{\text{BM}}(\mathbf{t}) = \begin{bmatrix} \text{zeros}(1, 21) \\ \mathbf{F}_{\text{BMx}}(\mathbf{t}) \\ \mathbf{F}_{\text{BMy}}(\mathbf{t}) \\ \text{zeros}(1, 1) \end{bmatrix} \quad (19)$$

and

$$\mathbf{F}_{\text{ext}}(\mathbf{t}) = \begin{bmatrix} \text{zeros}(1, 2) \\ \mathbf{C}_{\text{M}}(\mathbf{t}) \\ \text{zeros}(1, 18) \\ \mathbf{C}_{\text{R}}(\mathbf{t}) \end{bmatrix} \quad (20)$$

The vectors $\mathbf{F}_{\text{BMx}}(\mathbf{t})$ and $\mathbf{F}_{\text{BMy}}(\mathbf{t})$ are point vectors corresponding to the values of the forces exerted by the balls on the mill drum and calculated separately using Yade software. $[\mathbf{M}]_{(24 \times 24)}$ and $[\mathbf{K}]_{(24 \times 24)}$ are, respectively, the mass matrix and the stiffness matrix of the gearbox, obtained by assembling the elementary matrices of the system components. $\{\mathbf{X}\}_{(24 \times 1)}$ is the vector of degrees of freedom of the global system. $[\mathbf{C}]_{(24 \times 24)}$ represents the global damping matrix, calculated proportionally to the mass matrix $[\mathbf{M}]_{(24 \times 24)}$ and the average value of the stiffness matrix $[\mathbf{K}]_{(24 \times 24)}$ as follows:

$$\mathbf{C} = \alpha[\mathbf{M}] + \beta[\mathbf{K}] \quad (21)$$

$\alpha$ and $\beta$ are two real coefficients calculated so that the damping is viscous. According to [31], all resolution schemes for which $\alpha = 1/2$ and $\beta = 1/4$ are unconditionally stable. In this work, the resolution of the system of equations was conducted using the standard Newmark integration scheme. The mill torque (resistive torque) $\mathbf{C}_{\text{R}}$, the driving torque $\mathbf{C}_{\text{M}}$, and the speed of rotation $N_r$ were constants. Only the contact forces of the balls on the mill drum were variable. According to the parameters of the gear transmission systems summarized in Table 3, MATLAB software was used to solve Equation (14).

## 4. Results and Discussion

### 4.1. Effects of Grinding Forces

#### 4.1.1. Principle

In this section, the influence of the grinding forces is investigated. Three operating parameters were varied to understand the manifestations of the impulse forces in each case. The study was based on the observation of the load profile inside the drum, and interpretation of its influence by observing the amplitude of the signals in the vibrations time series and on the spectrogram. The maximum drop height of the balls $H_f$ was analyzed in all three cases. It was assumed that a high drop height correlated with a high impact force amplitude and consequently a high noise in the signal. A laboratory ball mill, whose characteristics were taken from the work of [4] and collected in Table 2, was used as a reference model. The simulations were performed according to the experimental design represented in Table 4, using Yade software. Yade (Yet another dynamic engine) is an open-source software tool based on the discrete element method (DEM), which uses object-oriented programming techniques [32]. The software program is written in C++ and has a
Python interface for high flexibility and performance. This makes the simulation of millions of particles interacting with each other possible. Yade runs on Linux operating systems and is available from the official repositories of common distributions such as Ubuntu. The source code is freely available on GitHub and is licensed under the GPL2 open-source license [22]. (Figure 7) depicts the results of a simulation of ball mill grindings forces by showing the evolution of contact forces between the balls and the liner, for three different levels of filling rate $F_r$. By analyzing the curve, we can observe that the average value of the applied forces increases with the filling rate.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Contact forces' time evolution, for three levels of filling rate: (a) $F_r$ = 15%, CAF = 496.5 N, (b) $F_r$ = 25%, CAF = 1772.5 N, (c) $F_r$ = 35% CAF = 5145.8 N, with $N_r$ = 0.3 $C_s$ and $D_b$ = 50 mm.</figcaption></figure>

The following subsections present the results of the effects of impulsions forces on the vibration behavior of the transmission for the three operating parameters.

### 4.1.2. Effect of the Filling Rate

Figure 8 shows the variation of the load profile inside the mill for different filling rates $F_r$ of 5%, 15%, 25%, and 35%. For a rotational speed $N_r$ of 30 rpm, and a ball diameter of 50 mm, the drop height of the balls $H_f$ increases with the increasing filling rate. For a filling rate of 15% (see Figure 8a), the maximum drop height is $H_f$ = 0.6262 m while for $F_r$ = 35%, $H_f$ = 0.8943 m. This increase in the drop height leads to an increase in the amplitude of the signals as can be seen in Figure 8b, in which the vibratory signature corresponding to $F_r$ = 35% presents a higher amplitude compared to that corresponding to $F_r$ = 15%.

The characteristics of the defect appear in this case with a severity of $\delta_f$ = 45%. The vibrations of the gearbox are affected in different ways by the increase or decrease of the filling rate. Table 6 presents a synthesis of the results for 16 similar observations, corresponding to four levels of crack severity, for four levels of filling rate. It provides information on the prevalence of different components in the vibration signal. For a filling rate $F_r$ of 5% and 15%, the mask noise $M_n$ prevails over the impulsive components of the defect IDC and over the impulsive noise $P_i$, when the defect degradation rate $\delta_f$ is between 0% and 30%. However, the impulsive noise is more dominant when the filling rate increases, i.e., for $F_r$ = 25% to $F_r$ = 35%.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Effects of filling rate: (a) load profile 1. $F_r$ = 5%, 2. $F_r$ = 15%, 3. $F_r$ = 25%, 4. $F_r$ = 35%, with $N_r$ = 0.6 $C_s$ and $D_b$ = 50 mm; (b) affected vibrations 1. time series of simulated signals $F_r$ = 35% and $F_r$ = 15%; 2. spectrogram of signal for $N_r$ = 0.6 $C_s$, $D_b$ = 50 mm.</figcaption></figure>

<table><caption>**Table 6.** Prevalence of the different phenomena in the signal due to the filling rate (Fr).</caption><thead><tr><th>Crack Level</th><th>F<sub>r</sub> = 5%</th><th>F<sub>r</sub> = 15%</th><th>F<sub>r</sub> = 25%</th><th>F<sub>r</sub> = 35%</th></tr></thead><tbody><tr><td>0%</td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>P<sub>i</sub></td><td>P<sub>i</sub></td></tr><tr><td>15%</td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>P<sub>i</sub></td><td>P<sub>i</sub></td></tr><tr><td>30%</td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>P<sub>i</sub></td><td>P<sub>i</sub></td></tr><tr><td>45%</td><td>M<sub>n</sub> + IDC</td><td>M<sub>n</sub> + IDC</td><td>M<sub>n</sub> + IDC</td><td>IDC + P<sub>i</sub></td></tr></tbody></table>

### 4.1.3. Effect of Rotation Speed (N)

The rotation speed shows two behaviors of the load profile inside the mill (Figure 9). For a filling rate of 30% and a ball diameter of 50 mm, the drop height increases together with increasing rotation speed. For a speed corresponding to 0.6 $C_s$, the maximum drop height is $H_f$ = 0.6366 m in Figure 9(a2). The balls rise to the specified drop height, then fall and roll over each other. The load has a cascade effect, and in practice, this effect favors attritional crushing. This can result in a high percentage of indirect contact between balls and walls, as most of the balls fall on top of each other. Beyond 0.6 $C_s$, the maximum drop height is $H_f$ = 0.7351 m (Figure 9(a4)). In this case, the balls rise to the height of drop and fall on the drum wall; a cataract effect is observed. In practice, this behavior favors impact crushing and thus the presence of a direct contact between ball and wall. In both cases, the vibration signals of the gearbox are amplified by the increase of the rotation speed (Figure 9b). However, this increase of the amplitudes remains weak compared to the augmentation caused by the increase of the filling rate. Figure 9b shows a weak amplification of the signal recorded at a rotation speed of 0.3 $C_s$, and is compared to the signal recorded at a rotation speed of 0.9 $C_s$ (critical speed). The increase in speed subjects the balls to a force approaching centrifugal force, thus the balls tend to stick to the inner wall of the mill, reducing the impact of the contact forces. Therefore, impulses due to defects are more visible at higher speeds since the impulses due to shocks are lower.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Effects of rotation speed: (a) load profile 1. N<sub>r</sub> = 0.3 C<sub>s</sub>, 2. N<sub>r</sub> = 0.6 C<sub>s</sub>, 3. N<sub>r</sub> = 0.7 C<sub>s</sub>, 4. N<sub>r</sub> = 0.9 C<sub>s</sub>, with F<sub>r</sub> = 30%, and D<sub>b</sub> = 50 mm; (b) affected vibrations 1. time series of simulated signals for N<sub>r</sub> = 0.9 C<sub>s</sub> (black), and N<sub>r</sub> = 0.3 C<sub>s</sub> (blue); 2. spectrogram of signal for F<sub>r</sub> = 30%, D<sub>b</sub> = 50 mm and δ<sub>f</sub> = 30%.</figcaption></figure>

Table 7 presents a synthesis of the results for 16 similar observations, corresponding to four levels of crack severity, for four levels of rotation speed. All the simulated signals show the prevalence of the mask noise $M_n$ and an easier manifestation of fault impulses for a degradation rate of $\delta_f$ = 30 and 45%.

<table><caption>**Table 7.** Prevalence of the different phenomena in the signal due to the rotation speed N<sub>r</sub>.</caption><thead><tr><th>Crack Level</th><th>N<sub>r</sub> = 0.3 C<sub>s</sub></th><th>N<sub>r</sub> = 0.6 C<sub>s</sub></th><th>N<sub>r</sub> = 0.7 C<sub>s</sub></th><th>N<sub>r</sub> = 0.9 C<sub>s</sub></th></tr></thead><tbody><tr><td>0%</td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td></tr><tr><td>15%</td><td>P<sub>i</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td></tr><tr><td>30%</td><td>M<sub>n</sub></td><td>P<sub>i</sub></td><td>IDC</td><td>IDC</td></tr><tr><td>45%</td><td>IDC</td><td>IDC</td><td>P<sub>i</sub></td><td>M<sub>n</sub></td></tr></tbody></table>

#### 4.1.4. Effect of Balls Size

The reduction in ball size leads to a reduction in the drop height creating extra free space in the mill. Considering the same rotational speed and the same filling rate (Figure 10), the balls can move freely and create a direct contact with the wall. The decrease in ball size promotes a cataract effect. For balls of 120 mm in diameter, the maximum head was H<sub>f</sub> = 0.5966 m, while it was H<sub>f</sub> = 0.6911 m for balls with a diameter of 60 mm. Despite the decrease in head, we registered a small increase in amplitude of the recorded signal for a diameter of 120 mm, compared to the one recorded for a diameter of 60 mm, as presented in Figure 10a. The spectrogram of Figure 10b was obtained for a signal corresponding to F<sub>r</sub> = 30%, N<sub>r</sub> = 0.6 C<sub>s</sub>, δ<sub>f</sub> = 0%. It has a large amount of mask noise M<sub>n</sub>, compared to impulsive noise P<sub>i</sub>. Table 6 presents a synthesis of the results for 16 similar observations, corresponding to four levels of crack severity, for four levels of balls size. The reduction of ball size favors the presence of mask noise. The low representativeness of the impulses due to the grinding is normal and due to the filling rate of 30%, a value for which there is a domination of the mask noise (see Table 8). A reduction in the ball size can only lead to a reduction in the filling rate, which in turn leads to the presence of mask noise. The defect only appears starting at a severity level of 30%.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Effect of ball size: (a) load profile 1. $D_b$ = 30 mm, 2. $D_b$ = 60 mm, 3. $D_b$ = 90 mm, 4. $D_b$ = 120 mm, for $F_r$ = 30% and $N_r$ = 0.3 $C_s$; (b) affected vibrations 1. time series of simulated signals for $D_b$ = 120 mm (black) and $D_b$ = 30 mm (blue); 2. spectrogram of signal, for $N_r$ = 0.6 $C_s$, $F_r$ = 30% and of $\delta_f$ = 0%.</figcaption></figure>

<table><caption>**Table 8.** Prevalence of the different phenomena in the signal due to the size of the balls $D_b$.</caption><thead><tr><th rowspan="2">Crack Level</th><th colspan="4">Ball Size D<sub>b</sub></th></tr><tr><th>D<sub>b</sub> = 30 mm</th><th>D<sub>b</sub> = 60 mm</th><th>D<sub>b</sub> = 90 mm</th><th>D<sub>b</sub> = 120 mm</th></tr></thead><tbody><tr><td>0%</td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td></tr><tr><td>15%</td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td></tr><tr><td>30%</td><td>IDC</td><td>IDC</td><td>IDC</td><td>IDC</td></tr><tr><td>45%</td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>M<sub>n</sub></td><td>IDC + P<sub>i</sub></td></tr></tbody></table>

## 4.2. Case without Grinding Efforts

In this section, the system was simulated with a tooth crack defect through increasing the rate of defect degradation, corresponding to $\delta_f$ = 30%, and eccentricity defect of 1 µm, but without grinding efforts. The gearbox rotated at a speed of 500 RPM, and the corresponding meshing frequency was $Gmf$ = 550 Hz. Figure 11 shows the vibrations simulated at node 4 in the Y-axis direction. From the spectrum and spectrogram, respectively in Figures 11b and 12, one can see the frequency characteristics related to the meshing and the crack defect. The model simulation reveals the presence of the meshing frequency and its odd-order harmonics (Figure 12). The crack defect appears in the time series (Figure 11a) at a period $T_{crack}$ = 0.87 s, which corresponds on the spectrum to the presence of the side bands spaced at a multiple frequency of the defect frequency $F_f$ = 1.1 Hz. As shown in Figure 12, the defect is manifested by the presence of side bands around the $Gmf$ frequency separated by the rotation frequency $f_r$.
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Simulation results without grinding efforts: (a) time series signals, (b) spectrogram of the signal with a crack defect $\delta_f = 30\%$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Spectrum showing manifestation of the eccentricity defect.</figcaption></figure>

### 4.3. Case with Grinding Efforts

This section shows the simulation results of the signals recorded at node 4, after the gearbox was loaded with grinding forces according to the different operating parameters studied previously. The analysis of Tables 5–7 reveals three important phenomena that can be observed on the time series, the spectrogram (Figure 13b), and the signal spectrum (Figure 14). For these three parameters, four levels of manifestations, for four different levels of severity of crack defect were retained. The idea was to see for each selected parameter level how the crack defect manifested itself in the presence of grinding forces. The vibratory signal in the presence of grinding forces is characterized by the presence of Gaussian noise *Mn*, as well as a noncyclic and non-Gaussian impulsive noise distribution
noted $P_i$. This noise is manifested by noncyclic peaks in the time series, by the presence of high-energy components in the spectrogram (Figure 13b), and by a random apparition of several lines around the gear mesh frequency in the spectrum (Figure 14). The distribution of impulsive noise in the signal is random and does not correspond to any fault-related feature. These two types of noise are present in all signals subjected to grinding forces. In the case of the presence of these noises, the cyclic impulsive component related to the crack defect IDC (impulsive defect component) is hidden. However, in the presence of grinding forces, this defect is masked by the impulses of the balls. In addition, in the case of a fault detection using these features, the presence of its lateral lines due to pulses can strongly influence the value of indicators calculated in both the time and frequency domains.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Simulation results with grindings efforts on node 5: (a) times series signal (b) spectrogram of signal with a crack defect $\delta_f = 30\%$ and $F_r = 30\%$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Spectrum of simulated signal with $\delta_f = 30\%$ and $F_r = 30\%$.</figcaption></figure>
## 5. Conclusions

This paper presented an original approach to study the influence of grinding forces on the vibratory behavior of a ball mill gearbox. By inserting a discrete element model of a ball mill into the dynamic model of a gear transmission, the paper highlighted for the first time three parameters intrinsic to the grinding process, i.e., the filling rate, the rotation speed, and the ball size, as the cause of the impulsive noise in the vibration signal. This is of great importance for transmission systems mounted in the mining industry operating in a noisy and impulsive environment. Through the short-term Fourier transform (STFT) of the envelope signal, the results showed that the increase of the filling rate leads to an increase of the impulsive noise. This phenomenon masked the defect signature, which only manifested itself at a high degradation rate. In this case, the defect was only visible in the spectrogram at a high severity level, i.e., 30–45% of the degradation rate. An increase of the rotation speed led to a decrease in the impulsive noise and increased the risk of the defect manifestation. In fact, when the rotation speed was close to the critical speed of the mill, the balls, due to centrifugal force, adhered to the mill wall and led to a decrease in impulses. In this case, the defect was visible in the spectrogram at a mean severity level, i.e., 15–30% of the degradation rate. The decrease in ball size had the opposite effect to that of the filling rate. It led to a manifestation of Gaussian noise. In this case, the defect was only visible in the spectrogram at a mean severity level, i.e., 30% of the degradation rate. Although all three parameters each had a significant contribution to impulsive noise, filling rate was the parameter that showed a large prominence over the other two.

The results presented in this paper were obtained by considering a constant load and speed. However, the dynamics of the load inside the mill implies a fluctuating velocity and load, which was not considered. These operating conditions will be considered variables in future work. In addition, the presence of impulsive shocks in the signal also contributes to masking the fault signal. Thus, signal-processing methods must be implemented in order to separate the impulsive shocks of the fault from those caused by the balls falling against the mill drum. The present work was limited to a numerical study of the observed phenomena. However, an experimental validation of the proposed model will be an important step towards the confirmation of the simulation results.

**Author Contributions:** G.N.K.: conceptualization, formal analysis, methodology, software, writing—original draft and visualization; X.C., L.R., A.B. and B.K.: funding acquisition, project administration, resources, supervision, validation, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** Within his PhD thesis, G.N.K. benefits from a French government scholarship, through the French embassy in DRC, managed by CAMPUS FRANCE under the grant number 113423W.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. Gearing Matrix

The modeling assumptions made here ensure that the pinions are comparable to non-deformable solids, except at certain points along the normal to the contact. The knowledge in the elementary reference frame (X₀, Y₀, Z₀) of the generalized displacements of the attachment nodes of the gearing element, allows the calculation of the displacements at the primitive meshing points. The vector of displacements and rotations of the node Nᵢ in the reference R₀ is:

$$ \left\{ \vec{\text{dep}}(N_i) \right\}_{R_0} = \begin{vmatrix} u_i & & \theta_{xi} \\ v_i & (24) & \theta_{yi} \\ w_i & & \theta_{zi} \end{vmatrix} \quad (A1) $$
It is then possible to calculate the displacements at the pitch point for each wheel:

$$ \{\vec{\text{dep}}(I \in S_i)\}_{R_o} = \{\vec{\text{dep}}(N_i)\}_{R_o} + \{\vec{\text{IN}}_i\} \wedge \{\vec{\text{rot}}(N_i)\}_{R_o} \quad (A2) $$

With:

$$ \{\vec{\text{IN}}_i\}_{R_o} = \begin{vmatrix} \mu_i \\ \lambda_i \\ -d_i \end{vmatrix} \quad (A3) $$

$d_i$ is the algebraic distance $N_iO_i$ measured along the axis $Z_o$. The values of the coefficients $\mu_i$ and $\lambda_i$ depend on the gearing configuration:

$$ \text{Configuration 1 :} \quad \begin{vmatrix} \mu_1 = -R_{b1} \tan \alpha \\ \mu_2 = +R_{b2} \tan \alpha \\ \lambda_1 = R_{b1} \\ \lambda_2 = -R_{b2} \end{vmatrix} \quad (A4) $$

and

$$ \text{Configuration 2 :} \quad \begin{vmatrix} \mu_1 = -R_{b1} \tan \alpha \\ \mu_2 = +R_{b2} \tan \alpha \\ \lambda_1 = -R_{b1} \\ \lambda_2 = R_{b2} \end{vmatrix} \quad (A5) $$

The displacements at the primitive meshing points are therefore expressed as:

$$ \{\vec{\text{dep}}(I \in S_i)\}_{R_o} = \begin{vmatrix} \mu_i + \lambda_i + d_i \theta_{yi} \\ v_i - d_i \theta_{xi} - \mu_i \theta_{zi} \\ w_i + \mu_i \theta_{yi} - \lambda_i \theta_{xi} \end{vmatrix} \quad (A6) $$

The expression of the contact normal in the reference $R_o$ depends on the mesh configuration and the value of the basic propeller angle $\beta$. The configuration chosen for the static and dynamic models is related to the direction of the helix of the wheel 1. If the helix is to the right of wheel 1, the helix angle is “negative”; it is positive in the case of a left helix. The outgoing normal to the contact profile of wheel 1 has the following expression in the reference frame $R_o$:

$$ \text{Configuration 1 :} \quad \{\vec{n}_1\}_{R_o} = \begin{vmatrix} \cos(\beta) \\ 0 \\ \sin(\beta) \end{vmatrix} \quad (A7) $$

and

$$ \text{Configuration 2 :} \quad \{\vec{n}_1\}_{R_o} = \begin{vmatrix} \cos(\beta) \\ 0 \\ -\sin(\beta) \end{vmatrix} \quad (A8) $$

The deformation at the teeth is then written as:

$$ \Sigma = \{\vec{\text{dep}}(I \in S_1)\} \cdot \{\vec{n}_1\} - \{\vec{\text{dep}}(I \in S_2)\} \{\vec{n}_1\} \quad (A9) $$

Let $\{X_{eng}\}$ be the vector of generalized displacements of the attachment nodes of the outer “cylindrical gear element”:

$$ \{X_{\text{gear}}\}_{R_o} = \{x_p, y_p, z_p, \theta_{xp}, \theta_{yp}, \theta_{zp}, x_g, y_g, z_{2g}, \theta_{xg}, \theta_{yg}, \theta_{zg}\}^T \quad (A10) $$

In this work, three degrees of freedom for each node were considered:

$$ \{X_{\text{gear}}\}_{R_o} = \{x_p, y_p, \theta_{zp}, x_g, y_g, \theta_{zg}\}^T \quad (A11) $$
The deformation at the teeth can then take the form of:

$$ \Sigma = \{\sigma\}^{\mathrm{T}}\{\mathbf{X}_{\mathrm{eng}}\}_{\mathrm{R}_{0}} \qquad (\mathrm{A}12) $$

With $\sigma$ the geometric vector of the transmission:

$$ \sigma = \{ c \ 0 \ c \ -\lambda_1 \ \mu_1 s + d_1 c \ \lambda_1 c \ -c \ 0 \ -s \ \lambda_2 s \ -d_2 c \ -\mu_2 s \ -\lambda_2 c \} \quad (\mathrm{A}13) $$

To consider the retained degrees of freedom, the geometric vector can be put in the following form:

$$ \sigma = \{ c \ 0 \ \lambda_1 c \ -c \ 0 \ -\lambda_2 c \} \qquad (\mathrm{A}14) $$

$$ c = \cos \beta \qquad (\mathrm{A}15) $$

$$ s = \sin \beta \qquad (\mathrm{A}16) $$

The meshing interface is modeled by a spring of stiffness $k_m$, placed along the contact normal. The potential energy dissipated in this element is expressed as:

$$ \begin{aligned} E_p &= \frac{1}{2} k_m \cdot \Sigma^2 = \frac{1}{2} k_m (\{\sigma\}^{\mathrm{T}} \cdot \{\mathbf{X}_{\mathrm{gear}}\}_{\mathrm{R}_o})^2 \\ &= \frac{1}{2} k_m (\{\mathbf{X}_{\mathrm{gear}}\}_{\mathrm{R}_o} \{\sigma\}^{\mathrm{T}} \cdot \{\sigma\} \{\mathbf{X}_{\mathrm{eng}}\}_{\mathrm{R}_o}^{\mathrm{T}}) = \frac{1}{2} k_m (\{\mathbf{X}_{\mathrm{gear}}\}_{\mathrm{R}_o} [\mathbf{K}_{\mathrm{gear}}] \{\mathbf{X}_{\mathrm{gear}}\}_{\mathrm{R}_o}^{\mathrm{T}}) \end{aligned} \qquad (\mathrm{A}17) $$

$$ [\mathbf{K}_{\mathrm{gear}}]_{\mathrm{R}_o} = k_m \{\sigma\}^{\mathrm{T}} \{\sigma\} \qquad (\mathrm{A}18) $$

(A18) is the equivalent stiffness matrix of the element “external cylindrical gear” expressed in the reference $R_0$

## Appendix B. Shafts Matrix

The input and output shafts of the gearbox were modeled by finite elements. They were each discretized into four elements, which were modeled by the Timoshenko beam with two nodes and six degrees of freedom. In this work, three degrees of freedom per node ($x_i, y_i, \theta_i$) were retained to simulate the bending and torsion effects of the shaft, i.e., two translations and one rotation. Gyroscopic and centrifugal effects were neglected. The displacement fields of the shaft element were approximated using the nodal displacements and shape functions as follows:

$$ \left\{ \begin{aligned} x(z,t) &= \{N_1(z)\}^t \{\boldsymbol{\chi}_x^e(t)\} \\ y(z,t) &= \{N_2(z)\}^t \{\boldsymbol{\chi}_y^e(t)\} \\ \theta_z(z,t) &= \{N_3(z)\}^t \{\boldsymbol{\chi}_{\theta_z}^e(t)\} \end{aligned} \right. \qquad (\mathrm{A}19) $$

where $N_1$ and $N_2$ are cubic polynomials and $N_3$ is a linear one. $\{\boldsymbol{\chi}_x^e(t)\}$ And $\{\boldsymbol{\chi}_y^e(t)\}$ are nodal displacements in the (XZ) and (YZ) planes while $\{\boldsymbol{\chi}_{\theta_z}^e(t)\}$ is the nodal displacements in the rotational directions. The elementary mass matrix $M_e$ and stiffness matrix $K_e$ were obtained by determining the kinetic $E_c$ and potential $E_p$ energy of the shaft element:

$$ E_c = \frac{1}{2} \rho A \int_0^L (\dot{x}^2 + \dot{y}^2) dz + \frac{1}{2} I_o \int_0^L \dot{\theta}_z^2 dz \qquad (\mathrm{A}20) $$

$$ E_p = \frac{1}{2} \int_0^L G \cdot J \cdot \left( \frac{\partial \theta_z}{\partial z} \right)^2 dz + \frac{1}{2} \int_0^L E \cdot I_A \cdot \left( \left( \frac{\partial^2 x}{\partial z^2} \right)^2 + \left( \frac{\partial^2 y}{\partial z^2} \right)^2 \right) dz \qquad (\mathrm{A}21) $$

In these two equations, $\rho$, $G$, $E$ represent, respectively, the density, the shear modulus, and the Young's modulus of the considered beam. $A$ and $L$ represent the cross section and length of the beam, $I_A$ is the mass inertia per unit length of the beam, along one of the axis
X or Y. $I_0$ is the polar inertia of the beam along the Z axis and J is the polar moment of inertia of the beam.

$$Me = \begin{bmatrix} \frac{156}{420}\rho AL + \frac{36}{30}I_x & 0 & 0 & \frac{54}{420}\rho AL - \frac{36}{30}I_x & 0 & 0 \\ 0 & \frac{156}{420}\rho AL + \frac{36}{30}I_x & 0 & 0 & \frac{54}{420}\rho AL - \frac{36}{30}I_x & 0 \\ 0 & 0 & \frac{I_0 L}{3} & 0 & 0 & \frac{I_0 L}{3} \\ \frac{54}{420}\rho AL - \frac{36}{30}I_x & 0 & 0 & \frac{156}{420}\rho AL + \frac{36}{30}I_x & 0 & 0 \\ 0 & \frac{156}{420}\rho AL + \frac{36}{30}I_x & 0 & 0 & \frac{156}{420}\rho AL + \frac{36}{30}I_x & 0 \\ 0 & 0 & \frac{I_0 L}{3} & 0 & 0 & \frac{I_0 L}{3} \end{bmatrix} \quad (A22)$$

$$Ke = \begin{bmatrix} \frac{12EI_A}{L^3} & 0 & 0 & -\frac{12EI_A}{L^3} & 0 & 0 \\ 0 & \frac{12EI_A}{L^3} & 0 & 0 & -\frac{12EI_A}{L^3} & 0 \\ 0 & 0 & \frac{GJ}{L} & 0 & 0 & -\frac{GJ}{L} \\ -\frac{12EI_A}{L^3} & 0 & 0 & \frac{12EI_A}{L^3} & 0 & 0 \\ 0 & -\frac{12EI_A}{L^3} & 0 & 0 & \frac{12EI_A}{L^3} & 0 \\ 0 & 0 & -\frac{GJ}{L} & 0 & 0 & \frac{GJ}{L} \end{bmatrix} \quad (A23)$$

## References

1. Jeswiet, J.; Szekeres, A. Energy Consumption in Mining Comminution. Procedia CIRP **2016**, *48*, 140–145. [CrossRef]

2. Górczyk, M.; Krot, P.; Zimroz, R.; Ogonowski, S. Increasing Energy Efficiency and Productivity of the Comminution Process in Tumbling Mills by Indirect Measurements of Internal Dynamics—An Overview. Energies **2020**, *13*, 6735. [CrossRef]

3. Molina Vicuna, C.; Venegas, P.J.; Valenzuela, M.A. Use of Virtual Sensors for the Analysis of Forces Exerted by the Load Inside a Tumbling Mill. Lyon, France, Juill. 2019. Available online: https://hal.archives-ouvertes.fr/hal-02190273 (accessed on 16 February 2021).

4. Boemer, D.; Ponthot, J.-P. A generic wear prediction procedure based on the discrete element method for ball mill liners in the cement industry. Miner. Eng. **2017**, *109*, 55–79. [CrossRef]

5. Kalala, J.T.; Breetzke, M.; Moys, M.H. Study of the influence of liner wear on the load behaviour of an industrial dry tumbling mill using the Discrete Element Method (DEM). Int. J. Miner. Process. **2008**, *86*, 33–39. [CrossRef]

6. Airikka, P. Simple Continuous-time Identification Method for PID Controlled Crushing Plant Processes. IFAC Proc. Vol. **2012**, *45*, 60–65. [CrossRef]

7. Mukepe Kahilu, M.; Kinnaert, M.; Kalenga Ngoy, P.; Moanda Ndeko, J.-M. MODELLING, IDENTIFICATION AND SIMULATION OF A DOUBLE CLOSED-LOOP INDUSTRIAL WET GRINDING CIRCUIT. IFAC Proc. Vol. **2012**, *45*, 72–77. [CrossRef]

8. Ruuska, J.; Lamberg, P.; Leiviskä, K. Flotation model based on floatability component approach—PGE minerals case. IFAC Proc. Vol. **2012**, *45*, 19–24. [CrossRef]

9. Chaari, F.; Bartelmus, W.; Zimroz, R.; Fakhfakh, T.; Haddar, M. Gearbox Vibration Signal Amplitude and Frequency Modulation. Shock Vib. **2012**, *19*, 635–652. [CrossRef]

10. Bartelmus, W.; Zimroz, R. Vibration condition monitoring of planetary gearbox under varying external load. Mech. Syst. Signal Process. **2009**, *23*, 246–257. [CrossRef]

11. Baydar, N.; Ball, A. DETECTION OF GEAR DETERIORATION UNDER VARYING LOAD CONDITIONS BY USING THE INSTANTANEOUS POWER SPECTRUM. Mech. Syst. Signal Process. **2000**, *14*, 907–921. [CrossRef]

12. Farhat, M.H.; Hentati, T.; Chiementin, X.; Bolaers, F.; Chaari, F.; Haddar, M. Numerical model of a single stage gearbox under variable regime. Mech. Based Des. Struct. Mach. **2020**, 1–28. [CrossRef]

13. Yu, G.; Shi, N. Gear fault signal modeling and detection based on alpha stable distribution. In Proceedings of the 2012 International Symposium on Instrumentation & Measurement, Sensor Network and Automation (IMSN), Sanya, China, 25–28 August 2012; IEEE: Piscataway, NJ, USA, 2012; pp. 471–474. [CrossRef]

14. Wylomanska, A.; Zimroz, R.; Janczura, J.; Obuchowski, J. Impulsive Noise Cancellation Method for Copper Ore Crusher Vibration Signals Enhancement. IEEE Trans. Ind. Electron. **2016**, *63*, 5612–5621. [CrossRef]

15. Schmidt, S.; Zimroz, R.; Chaari, F.; Stephan Heyns, P.; Haddar, M. A Simple Condition Monitoring Method for Gearboxes Operating in Impulsive Environments. Sensors **2020**, *20*, 2115. [CrossRef] [PubMed]

16. Wodecki, J. Local damage detection based on vibration data analysis in the presence of Gaussian and heavy-tailed impulsive noise. Measurement **2020**, *169*, 108400. [CrossRef]

17. Wodecki, J. Time-Varying Spectral Kurtosis: Generalization of Spectral Kurtosis for Local Damage Detection in Rotating Machines under Time-Varying Operating Conditions. Sensors **2021**, *21*, 3590. [CrossRef]

18. Hebda-Sobkowicz, J.; Zimroz, R.; Pitera, M.; Wylomańska, A. Informative frequency band selection in the presence of non-Gaussian noise—A novel approach based on the conditional variance statistic with application to bearing fault diagnosis. Mech. Syst. Signal Process. **2020**, *145*, 106971. [CrossRef]
19. Michalak, A.; Wodecki, J.; Drozda, M.; Wyłomańska, A.; Zimroz, R. Model of the Vibration Signal of the Vibrating Sieving Screen Suspension for Condition Monitoring Purposes. *Sensors* **2020**, *21*, 213. [CrossRef]

20. Wodecki, J.; Michalak, A.; Wyłomańska, A.; Zimroz, R. Influence of non-Gaussian noise on the effectiveness of cyclostationary analysis—Simulations and real data analysis. *Measurement* **2021**, *171*, 108814. [CrossRef]

21. Wodecki, J.; Kruczek, P.; Bartkowiak, A.; Zimroz, R.; Wyłomańska, A. Novel method of informative frequency band selection for vibration signal using Nonnegative Matrix Factorization of spectrogram matrix. *Mech. Syst. Signal Process.* **2019**, *130*, 585–596. [CrossRef]

22. Haustein, M.; Gladkyy, A.; Schwarze, R. Discrete element modeling of deformable particles in YADE. *SoftwareX* **2017**, *6*, 118–123. [CrossRef]

23. Overview—Yade, 3rd ed.; Documentation. Available online: http://yade-dem.org/doc/ (accessed on 7 March 2022).

24. Ashrafizadeh, H.; Ashrafizaadeh, M. Influence of processing parameters on grinding mechanism in planetary mill by employing discrete element method. *Adv. Powder Technol.* **2012**, *23*, 708–716. [CrossRef]

25. Hosseini, P.; Martins, S.; Martin, T.; Radziszewski, P.; Boyer, F.R. Acoustic emissions simulation of tumbling mills using charge dynamics. *Miner. Eng.* **2011**, *24*, 1440–1447. [CrossRef]

26. Gomez, J.L.; Khelf, I.; Bourdon, A.; André, H.; Rémond, D. Angular modeling of a rotating machine in non-stationary conditions: Application to monitoring bearing defects of wind turbines with instantaneous angular speed. *Mech. Mach. Theory* **2019**, *136*, 27–51. [CrossRef]

27. ISO No. 6336-1/2019; Calculation of Load Capacity of Spur and Helical Gears Part1: Basic Principles, Introduction and General Influence Factors 2019. International Organization for Standardization: Geneva, Switzerland, 2019. Available online: https://www.boutique.afnor.org/fr-fr/norme/iso-633612006/calcul-de-la-capacite-de-charge-des-engrenages-cylindriques-a-dentures-droi/xs018277/112945 (accessed on 10 March 2022).

28. Sghaier, E.; Bourdon, A.; Rémond, D.; Dion, J.-L.; Peyret, N. Coupled bending torsional vibrations of non-ideal energy source rotors under non-stationary operating conditions. *Int. J. Mech. Sci.* **2019**, *163*, 105155. [CrossRef]

29. Zhang, L.; Wu, C. Discrete element analysis of normal elastic impact of wet particles. *Powder Technol.* **2020**, *362*, 628–634. [CrossRef]

30. O'Sullivan, C.; Bray, J.D. Selecting a suitable time step for discrete element simulations that use the central difference time integration scheme. *Eng. Comput.* **2004**, *21*, 278–303. [CrossRef]

31. Parker, R.G.; Vijayakar, S.M.; Imajo, T. NON-LINEAR DYNAMIC RESPONSE OF A SPUR GEAR PAIR: MODELLING AND EXPERIMENTAL COMPARISONS. *J. Sound Vib.* **2000**, *237*, 435–455. [CrossRef]

32. Kozicki, J.; Donzé, F.V. YADE-OPEN DEM: An open-source software using a discrete element method to simulate granular material. *Eng. Comput.* **2009**, *26*, 786–805. [CrossRef]