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.