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 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 tenant compte :
- de la puissance mécanique à transmettre,
- de la résistance mécanique de l’acier utilisé,
- de la vitesse de rotation en sortie.
On doit donc déterminer :
- le rayon des roues,
- le nombre de dents,
- la largeur des dents,
- le coefficient de sécurité de l’ensemble .
Mais si les constantes physiques et les propriétés des matériaux sont en unités SI, les données du problème sont en unités impériales (dimensions en pouces) et en unités SI (puissance en watts). On a donc toutes les chances de se prendre les pieds dans le tapis en déroulant le calcul, et la NASA perd une sonde spatiale par décennie à cause d’erreurs de conversions d’unités de ce type.
On veut donc assurér la cohérence des résultats au cours du calcul en gardant la traçabilité des grandeurs et en évitant les erreurs de conversion.
1from pylab import *
2import quantities as u
3from tabulate import tabulate
4
5def inv(i):
6 return tan(i) - i
7
8### Données
9HB = 250. # dureté Rockwell B de l'acier
10p_d = 8 * 1/u.inch # pas diamétral de l'engrenage
11Delta_C = 0.1 * u.inch # tolérance max de l'entraxe
12R = 0.999 # fiabilité
13t = 1E5 # nombre de cycles de service (durée de vie utile)
14
15C = 5.5 * u.inch
16n_1 = 1200.0 * u.rpm # taux de rotation pignon (moteur)
17n_2 = 475.0 * u.rpm # taux de rotation roue (charge)
18P_u = 7.5 * u.kilowatt # puissance mécanique transmise par l'engrenage
19Q_v = 9. # facteur de fabrication
20S_p = 80 * u.megapascal # Contrainte plastique admissible de l'acier (non réversible)
21FS = 2. # Facteur de sécurité
22
23phi = (20. * u.deg).rescale('rad') # angle de contact
24
25### Calculs des rayons primitifs approchés
26i = n_1 / n_2
27r_1 = C / (i + 1)
28r_2 = i * r_1
29d_1 = 2. * r_1
30d_2 = 2. * r_2
31
32print """Calculs des rayons primitifs approchés
33----------------------
34 i = %s | Rapport de réduction
35 r_1 = %s | Rayon pignon
36 r_2 = %s | Rayon roue
37""" % (i, r_1, r_2)
38
39### Calcul du module
40W_t = 60 * P_u / (pi * d_1.rescale('m') * n_1)
41m = 2.34 * sqrt(W_t / (Q_v * S_p))
42
43print """Module
44----------------------
45 W_t = %s | chargement tangentiel
46 m = %s | module
47""" % (W_t.rescale('N'), m.rescale('mm') )
48
49### Calcul du nombre de dents
50Z_1 = p_d * d_1
51Z_2 = p_d * d_2
52
53print """Nombre approximatif de dents
54----------------------
55 Z_1 = %s | dents pignon
56 Z_2 = %s | dents roue
57""" % (Z_1, Z_2)
58
59Z_1_m = array([[floor(Z_1), ceil(Z_1)]])
60Z_2_m = array([[floor(Z_2)], [ceil(Z_2)]])
61i_m = (1. / Z_1_m) * Z_2_m
62C_d = 0.5 * Z_1_m / p_d + 0.5 * Z_2_m / p_d
63
64Z_1 = Z_1_m[1]
65Z_2 = Z_2_m[0,0]
66C = C_d[0,1]
67i = i_m[0,1]
68
69print """Sélection de la combinaison pignon/roue
70----------------------
71Rapports de réduction
72%s
73Entraxes
74%s
75Combinaison choisie
76 Z_1 = %i dents
77 Z_2 = %i dents
78 C = %s | entraxe
79 i = %s | rapport de transmission
80""" % (tabulate(column_stack((Z_2_m,i_m)), headers=Z_1_m.astype(int), tablefmt="fancy_grid"),
81 tabulate(column_stack((Z_2_m,C_d)), headers=Z_1_m.astype(int), tablefmt="fancy_grid"),
82 Z_1, Z_2, C, i
83 )
84
85### Dimension des dents
86a = 1. / p_d
87b = 1.25 / p_d
88h = a + b
89p_c = pi / p_d
90t = pi / (2. * p_d)
91r_f = 0.35 / p_d
92
93print """Dimension des dents
94----------------------
95 a = %s | Saillie
96 b = %s | Creux
97 h = %s | Hauteur
98 p_c = %s | Pas primitif
99 t = %s | Épaisseur au cercle primitif
100 r_f = %s | Rayon à la racine
101""" % (a, b , h, p_c, t, r_f)
102
103class engrenage:
104 def __init__(self, Z, p_d, a, b, phi):
105 self.d = Z / p_d
106 self.r = self.d / 2
107 self.d_a = self.d + 2*a
108 self.a = (self.d_a - self.d)/2
109 self.d_c = self.d - 2*b
110 self.r_b = self.r * cos(phi)
111 self.p_b = 2*pi*self.r_b /Z
112
113### Dimensions du pignon
114pignon = engrenage(Z_1, p_d, a, b, phi)
115
116print """Dimensions du pignon
117----------------------
118 d_1 = %s | Diamètre primitif
119 d_a_1 = %s | Diamètre saillie
120 d_c_1 = %s | Diamètre creux
121 r_b_1 = %s | Rayon de base
122 p_b = %s | Pas de base
123""" % (pignon.d, pignon.d_a, pignon.d_c, pignon.r_b, pignon.p_b)
124
125# Dimensions de la roue
126roue = engrenage(Z_2, p_d, a, b, phi)
127
128print """Dimensions de la roue
129----------------------
130 d_2 = %s
131 d_a_2 = %s
132 d_c_2 = %s
133 r_b_2 = %s
134 p_b = %s
135""" % (roue.d, roue.d_a, roue.d_c, roue.r_b, roue.p_b)
136
137### Vérification du contact et interférence
138m = ( 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
139g = (Z_1**2 + 2 * Z_1 * Z_2) * sin(phi)**2 - 4*Z_2 - 4
140print """Contact et interférences
141----------------------
142 m = %s
143 g(Z_1, Z_2) = %s
144""" % (m, g)
145
146### Paramètres d'opération
147C_dot = C + Delta_C
148r_dot_1 = C_dot / (1 + i)
149r_dot_2 = i * r_dot_1
150phi_dot = arccos(C * cos(phi) / C_dot)
151t_dot_1 = r_dot_1 * t / r_1 - 2 * r_dot_1 * (inv(phi_dot) - inv(phi))
152t_dot_2 = r_dot_2 * t / r_2 - 2 * r_dot_2 * (inv(phi_dot) - inv(phi))
153B = 2 * C_dot * (inv(phi_dot) - inv(phi))
154
155print """Paramètres d'opération
156----------------------
157 r_dot_1 = %s
158 r_dot_2 = %s
159 phi_dot = %s
160 t_dot_1 = %s
161 t_dot_2 = %s
162 B = %s
163""" % (r_dot_1,
164 r_dot_2,
165 phi_dot.rescale('deg'),
166 t_dot_1,
167 t_dot_2,
168 B)
169
170### Résistance en flexion
171S_fb_dot = (6235 + 174 * HB - 0.126*HB**2)*u.lb / u.inch**2
172K_L = 1.5
173K_T = 1
174K_R = 1.25
175S_fb = S_fb_dot * K_L / (K_T * K_R)
176sigma_b = S_fb / FS
177
178# Équation de Lewis
179W_T = 33000 * P_u.rescale('hp') / (pi * d_1.rescale('feet') * n_1) # lb
180
181J = 0.37
182K_v = 0.837
183K_m = 1.6
184F = W_T * K_m * p_d / (K_v * sigma_b * J)
185
186print """Résistance en flexion
187----------------------
188 S_fb = %s
189 sigma_b = %s
190 W_T = %s [lb]
191 F = %s [in]
192""" % (S_fb, sigma_b, W_T, F)
193
194### Largeur de dents
195C_L = 1.138
196C_H = 1
197C_T = 1
198C_R = 1.25
199S_dot_fc = (27000 + 364 * HB) * u.lb / u.inch**2
200
201S_fc = C_L * C_H * S_dot_fc / (C_T * C_R)
202
203C_p = 2300
204C_s = 1
205C_a = 1
206C_m = K_m
207C_v = K_v
208C_f = 1
209I = i * cos(phi) * sin(phi) / (2 * (i+1))
210sigma_c = S_fc / sqrt(FS)
211
212F_dot = C_p**2 * W_T/ (sigma_c**2 * I * d_1) * C_a * C_m * C_s * C_f / C_v
213
214print """Largeur des dents
215----------------------
216 S_dot_fc = %s
217 S_fc = %s
218 sigma_c = %s
219 F_dot = %s [in]
220""" % (S_dot_fc, S_fc, sigma_c, F_dot)</pre>Résultat
1Calculs des rayons primitifs approchés
2----------------------
3 i = 2.52631578947 dimensionless | Rapport de réduction
4 r_1 = 1.55970149254 in | Rayon pignon
5 r_2 = 3.94029850746 in | Rayon roue
6
7Module
8----------------------
9 W_t = 14386.2491373 N | chargement tangentiel
10 m = 10.4598004135 mm | module
11
12Nombre approximatif de dents
13----------------------
14 Z_1 = 24.9552238806 dimensionless | dents pignon
15 Z_2 = 63.0447761194 dimensionless | dents roue
16
17Sélection de la combinaison pignon/roue
18----------------------
19Rapports de réduction
20╒════╤═════════╤══════╕
21│ │ 24 │ 25 │
22╞════╪═════════╪══════╡
23│ 63 │ 2.625 │ 2.52 │
24├────┼─────────┼──────┤
25│ 64 │ 2.66667 │ 2.56 │
26╘════╧═════════╧══════╛
27Entraxes
28╒════╤════════╤════════╕
29│ │ 24 │ 25 │
30╞════╪════════╪════════╡
31│ 63 │ 5.4375 │ 5.5 │
32├────┼────────┼────────┤
33│ 64 │ 5.5 │ 5.5625 │
34╘════╧════════╧════════╛
35Combinaison choisie
36 Z_1 = 25 dents
37 Z_2 = 63 dents
38 C = 5.5 in | entraxe
39 i = 2.52 | rapport de transmission
40
41Dimension des dents
42----------------------
43 a = 0.125 in | Saillie
44 b = 0.15625 in | Creux
45 h = 0.28125 in | Hauteur
46 p_c = 0.392699081699 in | Pas primitif
47 t = 0.196349540849 in | Épaisseur au cercle primitif
48 r_f = 0.04375 in | Rayon à la racine
49
50Dimensions du pignon
51----------------------
52 d_1 = 3.125 in | Diamètre primitif
53 d_a_1 = 3.375 in | Diamètre saillie
54 d_c_1 = 2.8125 in | Diamètre creux
55 r_b_1 = 1.46826971998 in | Rayon de base
56 p_b = 0.369016429262 in | Pas de base
57
58Dimensions de la roue
59----------------------
60 d_2 = 7.875 in
61 d_a_2 = 8.125 in
62 d_c_2 = 7.5625 in
63 r_b_2 = 3.70003969434 in
64 p_b = 0.369016429262 in
65
66Contact et interférences
67----------------------
68 m = 1.70177405756 dimensionless
69 g(Z_1, Z_2) = 185.591113613 dimensionless
70
71Paramètres d'opération
72----------------------
73 r_dot_1 = 1.59090909091 in
74 r_dot_2 = 4.00909090909 in
75 phi_dot = 22.644361906 deg
76 t_dot_1 = 0.177860802115 in
77 t_dot_2 = 0.143285601241 in
78 B = 0.0789093806371 in
79
80Résistance en flexion
81----------------------
82 S_fb = 50232.0 lb/in**2
83 sigma_b = 25116.0 lb/in**2
84 W_T = 338.68017762 hp/(ft*rpm) [lb]
85 F = 0.557342411117 in*hp/(lb*ft*rpm) [in]
86
87Largeur des dents
88----------------------
89 S_dot_fc = 118000.0 lb/in**2
90 S_fc = 107427.2 lb/in**2
91 sigma_c = 75962.5016039 lb/in**2
92 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.
