Un des inconvénients de la programmation « mathématique » de grandeurs physiques est la difficulté de suivre les grandeurs utilisées au fil des variables et des calculs et donc d’assurer la cohérence du système d’unité et des multiples et sous-multiples utilisés. Habituellement, ceci passe par un lourd système de commentaires, peu efficace pour détecter les erreurs. Il apparaît donc pertinent de déclarer explicitement les unités des variables utilisées afin d’en assurer la traçabilité.
Cette pratique offre trois autres avantages :
- la possibilité de convertir à la volée les grandeurs d’un système d’unité à l’autre (vital en Amérique du Nord où le système métrique cohabite avec le système impérial et provoque un crash de sonde spatial par décennie à cause de soucis d’unités)
- la possibilité de débugguer un algorithme par les dimensions des résultats et les erreurs renvoyées en cas de tentatives d’opération illicites sur des grandeurs différentes
- la possibilité de programmer avec des grandeur sans se soucier de les convertir à chaque étape dans le bon multiple.
Il existe différents modules Pythons pour venir à bien de cette tâche, présentés dans la vidéo ci-dessus, chacun ayant leurs forces et leurs faiblesses, et une intégration plus ou moins fine avec Numpy, qui est la base du calcul numérique sous Python. J’ai choisi – comme l’auteur de la présentation – quantities pour ce faire, et je vais vous présenter un cas réel commenté : la conception d’un engrenage (dimension et résistance des matériaux)
Exemple et code source
L’exemple ci-dessous est typique de l’ingénierie nord-américaine : on doit dimensionner un train d’engrenages en fonction de puissance à transmettre et de la vitesse de sortie, puis calculer le coefficient de sécurité de l’ensemble (contrainte mécanique dans l’acier). On doit déterminer le rayon des roues, le nombre de dents, et la largeur des dents. Mais si les constantes physiques et les propriétés des matériaux sont en SI, les données du problème sont en impérial (dimensions en pouces) et en SI (puissance en watts). On veut donc gérer la cohérence des résultats en gardant la traçabilité des grandeurs et en évitant les erreurs de conversion.
# -*- coding: utf-8 -*-
"""
Created on Tue Nov 24 02:39:46 2015
@author: aurelien
"""
from pylab import *
import quantities as u
from tabulate import tabulate
def inv(i):
return tan(i) - i
### Données
HB = 250 # dureté Rockwell de l'acier
p_d = 8 * 1/u.inch # pas diamétral de l'engrenage
Delta_C = 0.1 * u.inch # tolérance max de l'entraxe
R = 0.999 # fiabilité
t = 1E5 # nombre de cycles de service
C = 5.5 * u.inch
n_1 = 1200.0 * u.rpm # taux de rotation pignon (moteur)
n_2 = 475.0 * u.rpm # taux de rotation roue (charge)
P_u = 7.5 * u.kilowatt # puissance transmise par l'engrenage
Q_v = 9 # facteur de fabrication
S_p = 80 * u.megapascal # Contrainte plastique admissible de l'acier
FS = 2 # Facteur de sécurité
phi = (20 * u.deg).rescale('rad') # angle de contact
### Calculs des rayons primitifs approchés
i = n_1/n_2
r_1 = C/(i+1)
r_2 = i * r_1
d_1 = 2*r_1
d_2 = 2*r_2
print """Calculs des rayons primitifs approchés
----------------------
i = %s | Rapport de réduction
r_1 = %s | Rayon pignon
r_2 = %s | Rayon roue
""" % (i, r_1, r_2)
### Calcul du module
W_t = 60 * P_u / (pi * d_1.rescale('m') * n_1)
m = 2.34 * sqrt(W_t/(Q_v * S_p))
print """Module
----------------------
W_t = %s | chargement tangentiel
m = %s | module
""" % (W_t.rescale('N'), m.rescale('mm') )
### Calcul du nombre de dents
Z_1 = p_d * d_1
Z_2 = p_d * d_2
print """Nombre approximatif de dents
----------------------
Z_1 = %s | dents pignon
Z_2 = %s | dents roue
""" % (Z_1, Z_2)
Z_1_m = array([[floor(Z_1), ceil(Z_1)]])
Z_2_m = array([[floor(Z_2)], [ceil(Z_2)]])
i_m = (1/Z_1_m)* Z_2_m
C_d = 0.5*Z_1_m/p_d + 0.5*Z_2_m/p_d
Z_1 = Z_1_m[1]
Z_2 = Z_2_m[0,0]
C = C_d[0,1]
i = i_m[0,1]
print """Sélection de la combinaison pignon/roue
----------------------
Rapports de réduction
%s
Entraxes
%s
Combinaison choisie
Z_1 = %i dents
Z_2 = %i dents
C = %s | entraxe
i = %s | rapport de transmission
""" % (tabulate(column_stack((Z_2_m,i_m)), headers=Z_1_m.astype(int), tablefmt="fancy_grid"),
tabulate(column_stack((Z_2_m,C_d)), headers=Z_1_m.astype(int), tablefmt="fancy_grid"),
Z_1, Z_2, C, i
)
### Dimension des dents
a = 1.0 /p_d
b = 1.25/p_d
h = a+b
p_c = pi/p_d
t = pi/(2.0*p_d)
r_f = 0.35/p_d
print """Dimension des dents
----------------------
a = %s | Saillie
b = %s | Creux
h = %s | Hauteur
p_c = %s | Pas primitif
t = %s | Épaisseur au cercle primitif
r_f = %s | Rayon à la racine
""" % (a, b , h, p_c, t, r_f)
class engrenage:
def __init__(self, Z, p_d, a, b, phi):
self.d = Z / p_d
self.r = self.d / 2
self.d_a = self.d + 2*a
self.a = (self.d_a - self.d)/2
self.d_c = self.d - 2*b
self.r_b = self.r * cos(phi)
self.p_b = 2*pi*self.r_b /Z
### Dimensions du pignon
pignon = engrenage(Z_1, p_d, a, b, phi)
print """Dimensions du pignon
----------------------
d_1 = %s | Diamètre primitif
d_a_1 = %s | Diamètre saillie
d_c_1 = %s | Diamètre creux
r_b_1 = %s | Rayon de base
p_b = %s | Pas de base
""" % (pignon.d, pignon.d_a, pignon.d_c, pignon.r_b, pignon.p_b)
# Dimensions de la roue
roue = engrenage(Z_2, p_d, a, b, phi)
print """Dimensions de la roue
----------------------
d_2 = %s
d_a_2 = %s
d_c_2 = %s
r_b_2 = %s
p_b = %s
""" % (roue.d, roue.d_a, roue.d_c, roue.r_b, roue.p_b)
### Vérification du contact et interférence
m = ( sqrt( (r_1 + a_1)**2 - (r_1 * cos(phi)) **2) + sqrt( (r_2 + a_2)**2 - (r_2 * cos(phi))**2) - C * sin(phi)) / p_b
g = (Z_1**2 + 2 * Z_1 * Z_2) * sin(phi)**2 - 4*Z_2 - 4
print """Contact et interférences
----------------------
m = %s
g(Z_1, Z_2) = %s
""" % (m, g)
### Paramètres d'opération
C_dot = C + Delta_C
r_dot_1 = C_dot / (1 + i)
r_dot_2 = i * r_dot_1
phi_dot = arccos(C * cos(phi) / C_dot)
t_dot_1 = r_dot_1 * t / r_1 - 2 * r_dot_1 * (inv(phi_dot) - inv(phi))
t_dot_2 = r_dot_2 * t / r_2 - 2 * r_dot_2 * (inv(phi_dot) - inv(phi))
B = 2 * C_dot * (inv(phi_dot) - inv(phi))
print """Paramètres d'opération
----------------------
r_dot_1 = %s
r_dot_2 = %s
phi_dot = %s
t_dot_1 = %s
t_dot_2 = %s
B = %s
""" % (r_dot_1,
r_dot_2,
phi_dot.rescale('deg'),
t_dot_1,
t_dot_2,
B)
### Résistance en flexion
S_fb_dot = (6235 + 174 * HB - 0.126*HB**2)*u.lb / u.inch**2
K_L = 1.5
K_T = 1
K_R = 1.25
S_fb = S_fb_dot * K_L / (K_T * K_R)
sigma_b = S_fb / FS
# Équation de Lewis
W_T = 33000 * P_u.rescale('hp') / (pi * d_1.rescale('feet') * n_1) # lb
J = 0.37
K_v = 0.837
K_m = 1.6
F = W_T * K_m * p_d / (K_v * sigma_b * J)
print """Résistance en flexion
----------------------
S_fb = %s
sigma_b = %s
W_T = %s [lb]
F = %s [in]
""" % (S_fb, sigma_b, W_T, F)
### Largeur de dents
C_L = 1.138
C_H = 1
C_T = 1
C_R = 1.25
S_dot_fc = (27000 + 364 * HB)*u.lb / u.inch**2
S_fc = C_L * C_H * S_dot_fc / (C_T * C_R)
C_p = 2300
C_s = 1
C_a = 1
C_m = K_m
C_v = K_v
C_f = 1
I = i * cos(phi) * sin(phi) / (2 * (i+1))
sigma_c = S_fc / sqrt(FS)
F_dot = C_p**2 * W_T/ (sigma_c**2 * I * d_1) * C_a * C_m * C_s * C_f / C_v
print """Largeur des dents
----------------------
S_dot_fc = %s
S_fc = %s
sigma_c = %s
F_dot = %s [in]
""" % (S_dot_fc, S_fc, sigma_c, F_dot)
Sortie
Calculs des rayons primitifs approchés
----------------------
i = 2.52631578947 dimensionless | Rapport de réduction
r_1 = 1.55970149254 in | Rayon pignon
r_2 = 3.94029850746 in | Rayon roue
Module
----------------------
W_t = 14386.2491373 N | chargement tangentiel
m = 10.4598004135 mm | module
Nombre approximatif de dents
----------------------
Z_1 = 24.9552238806 dimensionless | dents pignon
Z_2 = 63.0447761194 dimensionless | dents roue
Sélection de la combinaison pignon/roue
----------------------
Rapports de réduction
╒════╤═════════╤══════╕
│ │ 24 │ 25 │
╞════╪═════════╪══════╡
│ 63 │ 2.625 │ 2.52 │
├────┼─────────┼──────┤
│ 64 │ 2.66667 │ 2.56 │
╘════╧═════════╧══════╛
Entraxes
╒════╤════════╤════════╕
│ │ 24 │ 25 │
╞════╪════════╪════════╡
│ 63 │ 5.4375 │ 5.5 │
├────┼────────┼────────┤
│ 64 │ 5.5 │ 5.5625 │
╘════╧════════╧════════╛
Combinaison choisie
Z_1 = 25 dents
Z_2 = 63 dents
C = 5.5 in | entraxe
i = 2.52 | rapport de transmission
Dimension des dents
----------------------
a = 0.125 in | Saillie
b = 0.15625 in | Creux
h = 0.28125 in | Hauteur
p_c = 0.392699081699 in | Pas primitif
t = 0.196349540849 in | Épaisseur au cercle primitif
r_f = 0.04375 in | Rayon à la racine
Dimensions du pignon
----------------------
d_1 = 3.125 in | Diamètre primitif
d_a_1 = 3.375 in | Diamètre saillie
d_c_1 = 2.8125 in | Diamètre creux
r_b_1 = 1.46826971998 in | Rayon de base
p_b = 0.369016429262 in | Pas de base
Dimensions de la roue
----------------------
d_2 = 7.875 in
d_a_2 = 8.125 in
d_c_2 = 7.5625 in
r_b_2 = 3.70003969434 in
p_b = 0.369016429262 in
Contact et interférences
----------------------
m = 1.70177405756 dimensionless
g(Z_1, Z_2) = 185.591113613 dimensionless
Paramètres d'opération
----------------------
r_dot_1 = 1.59090909091 in
r_dot_2 = 4.00909090909 in
phi_dot = 22.644361906 deg
t_dot_1 = 0.177860802115 in
t_dot_2 = 0.143285601241 in
B = 0.0789093806371 in
Résistance en flexion
----------------------
S_fb = 50232.0 lb/in**2
sigma_b = 25116.0 lb/in**2
W_T = 338.68017762 hp/(ft*rpm) [lb]
F = 0.557342411117 in*hp/(lb*ft*rpm) [in]
Largeur des dents
----------------------
S_dot_fc = 118000.0 lb/in**2
S_fc = 107427.2 lb/in**2
sigma_c = 75962.5016039 lb/in**2
F_dot = 1.65388227338 in**3*hp/(lb**2*ft*rpm) [in]
Commentaires
Sur la procédure de calcul, on remarque aux sections Dimension du pignon et Dimension de la roue que les calculs étant les mêmes, le code a été factorisé en passant par une classe. C’est l’une des beautés de la programmation scientifique orientée objet qui est grandement facilitée avec Python, et plus complexe – quoique possible – sous Matlab par exemple.
Sur les unités, on peut vérifier la cohérences des résultats, y compris sur les rapports de transmission sans dimension.
Cependant, on remaque un problème de calcul sur la variable W_T (section Résistance en flexion) dont la dimension est hp/(ft*rpm) en raison du facteur de correction 33000 sans unité précisée. Ainsi W_T est homogène à des livres force mais le passage dans cette unité via le module quantities introduit des constantes de conversion déjà prises en compte dans le 33 000. Le problème se propage à F. Note : La valeur 33000 a été reprise telle quelle d’un corrigé qui ne mentionnait pas sa dimension ni son origine (Car les formules de dimensionnement d’éléments de machines relèvent parfois de la sorcellerie).
À cette exception près (qui mériterait seulement un retour à la documentation), on remarque que le passage d’un système d’unités à l’autre est presque indolore.
