5.2 Conception d'un
modèle théorique d'estimation du temps de survie
Dans le chapitre sur l'analyse vibratoire nous avons vu que pour
surveiller une machine il fallait suivre l'évolution des indicateurs de
défauts de cette dernière. L'indicateur que nous avons
utilisé pour détecter l'apparition des premiers défauts du
roulement est la moyenne des cinq plus grandes valeurs absolues des
accélérations mesurées dans chaque fichier d'observation.
La moyenne a été utilisée pour réduire
l'effet du bruit. La formule qui définit cet indicateur est la
suivante :
Où accdésigne
l'accélération et i est le rang de la
donnée après avoir trié dans l'ordre décroissant
les valeurs absolues de l'accélération.
Cet indicateur doit être évalué dans
chaque instant de prélèvement, puis on doit générer
un graphe pour suivre son évolution. Pour atteindre cet objectif, nous
avons implémenté une fonction que nous avons baptisée
moyenne5.m dont l'algorithme est le suivant :
Début algorithme
Récupération du nombre de fichiers (instants
de prélèvement) du roulement traité
Pour i allant de 1 au nombre de fichiers
traités
Construction du nom du fichier à exploiter
Extraction des données contenues dans le fichier
traité
Récupération de la colonne
d'accélération horizontale
Récupération de la colonne
d'accélération verticale
Calcul des valeurs absolues des signaux temporels
Tri des valeurs absolues dans l'ordre
décroissant
Récupération des cinq premières valeurs
des tableaux
Calcul de la moyenne de ces valeurs
Insertion du résultat dans un tableau
Fin Pour
Sauvegarde des spectres des signaux dans des fichiers
EXCEL
Construction de l'axe des abscisses
Construction de la courbe d'évolution de
l'indicateur des deux signaux
Finalgorithme
Le code source complet de cette fonction sera fourni dans
l'annexe B. Nous allons revenir sur certaines spécificités de cet
algorithme vu qu'on retrouve quelques instructions qui apparaissent
déjà dans l'algorithme précédent.
Le calcul de la valeur absolue se fait à l'aide de la
fonction abs déjà contenue dans la
bibliothèque mathématique de MATLAB. Pour l'utiliser, il suffit
de lui passer en paramètre le vecteur contenant les nombres dont on
souhaite obtenir les valeurs absolues.
Le tri d'un tableau est réalisé à l'aide
de la fonction sort de MATLAB. On lui envoie un
vecteur en paramètre (dans notre cas le vecteur comporte les valeurs
absolues des signaux temporels d'accélération) et en sortie elle
nous renvoie le vecteur trié dans l'ordre croissant. Pour qu'elle puisse
trier le tableau dans l'ordre décroissant comme nous souhaitons le
faire, il faut lui ajouter un second paramètre, notamment la chaine de
caractères « descend ».
Lorsqu'on a un vecteur X dans MATLAB et qu'on
veut obtenir les n premières valeurs de ce
dernier, on utilise la syntaxe : X(1 : n). Dans
notre cas d'usage n = 5, car on souhaite
récupérer les cinq premières valeurs du vecteur
trié dans l'ordre décroissant. Ces cinq valeurs sont en effet les
plus grandes du vecteur. Pour finir, on calcul la moyenne de ces cinq valeurs
en utilisant la fonction
meanimplémentée dans MATLAB.
Pour générer le graphe d'évolution de
l'indicateur, nous avons utilisé la fonction plot
de MATLAB. Cette fonction permet de construire un graphe à
deux dimensions lorsqu'on lui fournit en paramètre l'axe des abscisses
dans un premier vecteur et l'axe des ordonnées dans un second vecteur.
Dans notre cas, en abscisse nous avions les numéros de chaque fichier
(instants de prélèvement) et en ordonnées nous avions les
valeurs de l'indicateur correspondant à chaque fichier d'observation.
Nous avons exécuté cet algorithme sur le
roulement bearing1_1 et nous avons obtenu le graphe de la figure
suivante :
Figure 35: indicateur appliqué aux
accélérations horizontales
Figure 36: indicateur appliqué aux
accélérations verticales
Le résultat obtenu sur le signal
d'accélération verticale n'est pas facile à exploiter. Par
contre celui obtenu sur le signal d'accélération horizontale
montre d'abord une évolution quasi constante de l'indicateur, et
à partir de l'observation #1218 on remarque une croissance graduelle. A
cet instant, le roulement entre dans la première phase (ou stade) de
dégradation.
A l'observation #2747 on observe une croissance brusque de
l'indicateur. Le roulement entre dans la deuxième de
dégradation.
La tendance de dégradation du roulement bearing1_1
est la même que celle observée dans le roulement
bearing1_3. Ce sont les deux seuls cas dans lesquels nous avons
observéune croissance progressive de l'indicateur après
l'apparition du premier défaut. Dans les autres roulements, l'indicateur
présente des successions de croissances et de décroissances
après l'apparition du premier défaut.
Nous allons donc concevoir le modèle d'estimation du
temps de survie à l'aide du roulement bearing1_1 et nous allons
l'appliquer sur le roulement bearing1_3.
Pour mettre en place le modèle de prédiction,
une approche à deux étapes a été adoptée.
La première étape consiste à estimer le
temps au bout duquel le roulement entre dans la deuxième phase
d'anomalie et la deuxième étape est d'estimer la durée
entre le début de la deuxième phase d'anomalie et l'instant
où le roulement est hors d'usage.
A partir du résultat obtenu à la figure 35,
l'idée est de générer une courbe d'ajustement
exponentielle à partir de l'instant de détection du
premier défaut (phase 1 d'anomalie) jusqu'au début de la phase 2
d'anomalie. Pour atteindre cet objectif, nous avons développé une
fonction dans MATLAB que nous avons
baptiséeexp_fit.m. Cette fonction obéit
à l'algorithme suivant :
Début algorithme
Récupération du nombre de fichiers (instants
de prélèvement) du roulement traité
Pour i allant de 1 au nombre de fichiers
traités
Construction du nom du fichier à exploiter
Extraction des données contenues dans le fichier
traité
Récupération de la colonne
d'accélération horizontale
Calcul des valeurs absolues des signaux temporels
Trie des valeurs absolues dans l'ordre
décroissant
Récupération des cinq premières valeurs
du tableau
Calcul de la moyenne de ces valeurs
Insertion du résultat dans un tableau
Fin Pour
Sauvegarde des spectres des signaux dans des fichiers
EXCEL
Troncature du tableau des valeurs de l'indicateur
Génération de la courbe d'ajustement
exponentiel
Construction de la courbe d'évolution de
l'indicateur et la courbe d'ajustement
Finalgorithme
Le code source complet de ce script sera fourni dans l'annexe
C de ce document. Revenons toutefois sur les étapes clés de cet
algorithme.
La troncature permet de récupérer les valeurs de
l'indicateur qui sont comprises entre la phase 1 et la phase 2 d'anomalie. Dans
MATLAB, lorsqu'on a un vecteur X et qu'on souhaite
récupérer les valeurs comprises entre les positions n
et m du vecteur, on utilise la syntaxe
suivante : X(n : m). Dans notre cas, X correspond au
tableau d'évolution des valeurs de l'indicateur appliqué à
l'accélération horizontale, n vaut 1218 et m vaut 2747.
Pour générer la courbe d'ajustement exponentiel
des données, on utilise la fonction fit
contenue dans MATLAB. Cette fonction prend trois
paramètres, à savoirl'axe des abscisses et des ordonnées
ainsi que la chaine de caractère
« exp1 » qui spécifie le
type d'ajustement qu'on applique aux données. Dans la dernière
version de MATLAB, un outil spécial a été
développé spécifiquement pour résoudre les
problèmes d'ajustement linéaire, logarithmique, exponentiel, etc.
Cet outil permet une manipulation beaucoup plus souple des courbes
d'ajustement.
Après l'exécution de ce script, nous avons
obtenu le graphe suivant :
Figure 37: courbe d'ajustement exponentiel appliquée
au roulement bearing1_1
Cette fonction renvoie également l'expression
littérale de la courbe d'ajustement exponentiel. Son équation est
de la forme :
f(x) = aebx
Le logiciel MATLAB nous indique que les coefficients
a et b sont estimés
dans un intervalle de confiance à 95% et valent respectivement
0,3757 et 0,00111. Au
final, l'expression littérale de l'équation de la courbe
d'ajustement exponentiel est :
f(x) = 0,3757e 0,00111x
L'image du début de la phase 1 d'anomalie par la
fonction d'ajustement est 1,453 comme on peut si bien
le lire sur le graphe de la figure 37, et celle du début de la phase 2
d'anomalie est 7,95. Leur rapport vaut :
On détermine ensuite le temps au bout duquel le
roulement s'abîme complètement après être
entré dans la phase 2 d'anomalie. Sachant que chaque instant de
prélèvement est équivalent à une durée de 10
secondes, ce temps vaut :
?t1 = tf - ti =
28030s - 27470s = 560s
On évalue également la durée entre les
phases 1 et 2 d'anomalie :
?t2 = tf - ti =
27470s - 12180s = 15290s
Les durées que nous venons de calculer correspondent
respectivement à celle de la première et de la deuxième
anomalie. Pour finir, on évalue leur rapport qu'on va appeler le
ratio d'anomalies :
| 
