Intégration numérique : les caractéristiques de section d'un polygone quelconque, à la manière des logiciels

Avant qu'une poutre soit vérifiée en flexion, en flèche ou en flambement, sa section doit être réduite à une poignée de nombres. Ces nombres sont les caractéristiques de section :

• L'aire A, qui reprend l'effort normal et le poids propre.
• Le centre de gravité (C_x, C_y), le point par rapport auquel toutes les autres caractéristiques sont mesurées.
• Le moment quadratique I_x, I_y (le moment d'inertie), qui fixe la rigidité de flexion et, via la charge d'Euler, la capacité au flambement.
• Le produit d'inertie I_xy, qui indique si x et y sont les axes principaux ou si la section va fléchir latéralement sous une charge verticale.
• Le module de flexion S = I / c et le rayon de giration r = √(I/A), les nombres dérivés que les vérifications de la norme utilisent directement.

Pour un rectangle ou un rond, vous trouvez ces valeurs dans une formule. Mais une section réelle, un PRS soudé, une lèvre formée à froid, un gousset à coin coupé, une cellule de pont, n'a pas de formule. Ce n'est qu'un contour. La question à laquelle répond cet article est celle que tout programme de structure doit résoudre en interne : à partir des seules coordonnées des coins d'une forme quelconque, comment obtenir A, le centre de gravité et l'inertie, exactement et vite ?

La réponse honnête est l'intégration numérique, mais pas celle que la plupart des gens imaginent. Le logiciel ne découpe pas l'intérieur en grille pour additionner de petits carrés. Il fait quelque chose de bien plus élégant, qui sort tout droit de l'analyse.

Ouvrez le cours d'analyse au chapitre des intégrales doubles et vous avez déjà les définitions. Sur la région R qu'occupe la section :

A = ∫∫_R dA

Q_x = ∫∫_R y dA, Q_y = ∫∫_R x dA (les moments statiques ; le centre de gravité vaut C_y = Q_x/A, C_x = Q_y/A)

I_x = ∫∫_R y² dA, I_y = ∫∫_R x² dA, I_xy = ∫∫_R xy dA

Ces expressions sont exactes et tout à fait générales, et si le lien entre une intégrale et une grandeur physique paraît abstrait, le cas de l'effort tranchant et du moment le rend concret : le moment est littéralement l'intégrale définie du diagramme d'effort tranchant. La difficulté ici est purement pratique : R est une région quelconque, il n'y a donc pas de primitive dans laquelle injecter des bornes. Un étudiant de première année poserait les bornes à la main pour une forme précise. Un programme ne le peut pas, car il ne connaît pas la forme avant l'exécution.

Il existe trois issues, et il vaut la peine de les nommer car on les confond souvent :

1. Maillage / Monte-Carlo : remplir l'intérieur de cellules ou de points aléatoires et sommer. Simple, marche toujours, mais approximatif et lent, et il faut l'intérieur, pas seulement le contour. C'est ce que fait un logiciel aux éléments finis pour la contrainte, pas ce qu'il fait pour les caractéristiques de section.
2. Quadrature de Gauss : échantillonner l'intégrande en quelques points astucieux avec des poids astucieux. Exacte pour les polynômes jusqu'à un certain degré, et comme y² est un polynôme de bas degré, une poignée de points suffit pour une forme à arêtes droites. C'est le cheval de trait des intégrales de rigidité d'élément d'un moteur aux éléments finis.
3. Théorème de Green : convertir l'intégrale d'aire en intégrale de contour, puis, pour un polygone, l'évaluer exactement comme une somme finie. Pas de maillage, pas d'échantillonnage, pas d'erreur. C'est ce qu'utilise un moteur de caractéristiques de section.

L'option 3 est celle qui transforme un problème non borné en une boucle d'une ligne, c'est donc celle que nous suivons.

Le théorème de Green est le résultat d'analyse qui relie une intégrale double sur une région à une intégrale de contour le long de sa frontière ∂R :

∫∫_R (∂M/∂x − ∂L/∂y) dA = ∮_∂R (L dx + M dy)

L'astuce consiste à choisir L et M pour que le membre de gauche devienne la caractéristique voulue. Pour l'aire, prenez L = 0, M = x, et le théorème donne A = ∮ x dy. Pour le moment quadratique par rapport à l'axe x, prenez le couple qui rend ∂M/∂x − ∂L/∂y = y², et I_x devient une intégrale de contour de termes en y². Le même geste marche pour chaque caractéristique de la liste.

Pourquoi est-ce un tel gain ? Parce que le contour d'un polygone est une suite de segments droits, et une intégrale de contour sur un segment droit se calcule en forme fermée une fois, symboliquement, pour toujours. Vous faites le calcul une seule fois sur une arête générique de (xᵢ, yᵢ) à (xᵢ₊₁, yᵢ₊₁), et à partir de là l'ordinateur n'additionne plus que l'algèbre obtenue. L'intégrale bidimensionnelle sur une région inconnue s'est réduite à une somme unidimensionnelle sur n arêtes.

Géométriquement, chaque arête balaie un triangle signé vers l'origine. Les arêtes du côté opposé à l'origine balaient une aire positive, celles du côté proche une aire négative, et les recouvrements s'annulent exactement, laissant la vraie région enfermée. Cette image du balayage signé est toute l'idée de la formule de la section suivante.

Menez l'intégrale de contour pour un polygone de sommets (xᵢ, yᵢ), listés dans l'ordre et fermé (le dernier sommet se reconnecte au premier), et chaque caractéristique se réduit à une somme d'une quantité par arête, le terme croisé :

crossᵢ = xᵢ · yᵢ₊₁ − xᵢ₊₁ · yᵢ

C'est exactement le numérateur d'un déterminant 2×2, l'aire signée du triangle formé par l'origine et l'arête i. L'aire vaut alors :

A = ½ Σ crossᵢ

C'est la formule du lacet (shoelace) (formule de l'aire de Gauss), nommée ainsi parce que le motif de multiplication croisée ressemble au laçage d'une chaussure. Elle est signée : listez les sommets dans le sens trigonométrique et A sort positif, dans le sens horaire et il sort négatif. Le logiciel exploite ce signe délibérément, positif pour un contour plein, négatif pour un trou, si bien qu'une section avec une cavité se traite en listant les sommets du trou dans l'autre sens et en additionnant.

Exemple résolu : un rectangle 200 × 300 mm

Plaçons les coins en (0,0), (200,0), (200,300), (0,300), dans le sens trigonométrique. Les quatre termes croisés :

• arête (0,0)→(200,0) : 0·0 − 200·0 = 0
• arête (200,0)→(200,300) : 200·300 − 200·0 = 60 000
• arête (200,300)→(0,300) : 200·300 − 0·300 = 60 000
• arête (0,300)→(0,0) : 0·0 − 0·300 = 0

A = ½ (0 + 60 000 + 60 000 + 0) = 60 000 mm², soit 200 × 300 au chiffre près. Deux des quatre arêtes n'ont rien apporté, car une arête située sur une droite passant par l'origine balaie un triangle nul. La formule n'en a cure ; elle additionne, tout simplement.

> Note CalcSteel : le moteur de section qui alimente les pages de profils exécute exactement cette boucle. Le terme croisé signé est la première ligne de la fonction ; tout le reste n'est que la même algèbre déclinée.

Le centre de gravité utilise le même terme croisé, pondéré par la coordonnée qui varie :

C_x = 1/(6A) Σ (xᵢ + xᵢ₊₁) · crossᵢ

C_y = 1/(6A) Σ (yᵢ + yᵢ₊₁) · crossᵢ

Le centre de gravité compte parce que le moment quadratique que veut la norme est mesuré par rapport au centre de gravité, pas par rapport à l'origine de votre dessin. Trompez-vous sur le centre de gravité et toute inertie qui en découle est fausse du terme de transport des axes parallèles.

Exemple résolu : un triangle rectangle, base 120 × hauteur 180

Sommets (0,0), (120,0), (0,180). L'aire par le lacet vaut ½ · 120 · 180 = 10 800 mm². En passant les mêmes sommets dans les sommes du centre de gravité, on obtient :

• C_x = 40 mm = base / 3
• C_y = 60 mm = hauteur / 3

ce qui est le centre de gravité du manuel pour un triangle rectangle, au tiers de chaque côté depuis le coin de l'angle droit. Remarquez que l'origine était posée sur un coin, pas sur le centre de gravité : la formule a trouvé le centre de gravité pour nous, et c'est tout l'intérêt. Vous n'avez jamais à savoir où se trouve le centre de gravité à l'avance.

Les moments quadratiques sortent de la même intégrale de contour, une étape algébrique de plus. Par rapport à l'origine du dessin :

I_x0 = 1/12 Σ (yᵢ² + yᵢ·yᵢ₊₁ + yᵢ₊₁²) · crossᵢ

I_y0 = 1/12 Σ (xᵢ² + xᵢ·xᵢ₊₁ + xᵢ₊₁²) · crossᵢ

I_xy0 = 1/24 Σ (xᵢ·yᵢ₊₁ + 2xᵢ·yᵢ + 2xᵢ₊₁·yᵢ₊₁ + xᵢ₊₁·yᵢ) · crossᵢ

Ces valeurs sont par rapport à l'origine, ce qui n'est presque jamais où vous les voulez. La dernière étape est le théorème des axes parallèles (le transport de Huygens/Steiner), utilisé à l'envers pour glisser l'axe de l'origine vers le centre de gravité :

I_x = I_x0 − A · C_y², I_y = I_y0 − A · C_x², I_xy = I_xy0 − A · C_x · C_y

Comme on soustrait la contribution du centre de gravité une seule fois, à la fin, tout le calcul reste un unique passage sur les arêtes suivi de trois soustractions. Rien n'est maillé, rien n'est échantillonné.

Retour au rectangle

Pour le rectangle 200 × 300, les sommes à l'origine donnent I_x0, et après le transport de C_y = 150 mm le résultat est I_x = 450 × 10⁶ mm⁴, exactement b·h³/12 = 200·300³/12. De même I_y = 200 × 10⁶ mm⁴. De là tombent les nombres dérivés qu'utilise la norme : S_x = I_x/c = 450×10⁶/150 = 3,0 × 10⁶ mm³, et r_x = √(I_x/A) = 86,6 mm. La voie polygonale reproduit la forme fermée jusqu'au dernier chiffre, et c'est bien le point : ce n'est pas une approximation de la formule, c'est la formule, dérivée une fois pour une arête générique.

Le produit d'inertie I_xy est la caractéristique dont personne ne se souvient dans la table de formules, car pour toute section admettant un axe de symétrie il est nul. Mais une cornière, un Z, un profil en U monosymétrique chargé hors axe, ou toute forme rognée ou dissymétrique a un I_xy non nul, et cela change la physique : une charge verticale produit une flexion latérale, car les axes fort et faible sont tournés hors de l'horizontale et de la verticale.

Le moteur obtient I_xy de la même somme en un passage montrée plus haut, puis trouve les axes principaux avec les relations du cercle de Mohr des inerties, qui ne sont que de la trigonométrie de niveau analyse, la rotation d'une forme quadratique :

R = √( ((I_x − I_y)/2)² + I_xy² ), I₁ = (I_x+I_y)/2 + R, I₂ = (I_x+I_y)/2 − R

θ = ½ · atan2(2 I_xy, I_y − I_x) situe l'axe.

Exemple résolu : une cornière à ailes égales L100 × 100 × 10

Tracez le L en six sommets, coin à l'origine, les deux ailes de 10 mm d'épaisseur. Un passage renvoie :

• A = 1 900 mm², centre de gravité en (28,68, 28,68), sur la ligne à 45° comme l'exige la symétrie.
• I_x = I_y = 1,80 × 10⁶ mm⁴ (égaux, parce que les ailes sont égales).
• I_xy = −1,066 × 10⁶ mm⁴, et le fait qu'il ne soit pas nul est toute l'histoire.

Les valeurs principales sont I₁ = 2,87 × 10⁶ mm⁴ et I₂ = 0,734 × 10⁶ mm⁴, à θ = 45°. Ainsi le vrai axe faible d'une cornière à ailes égales ne reprend qu'environ le quart de la rigidité de son axe fort, et les deux sont inclinés de 45° par rapport aux ailes. Un ingénieur qui vérifierait la flexion autour de l'axe de l'aile horizontale, avec I_x = 1,80 × 10⁶, se tromperait à la fois en grandeur et en direction. Le moteur polygonal vous livre I₂ et θ sans que rien de tout cela ne soit évident à partir du contour.

Voici l'exemple qui montre pourquoi cette méthode mérite confiance pour de vrais profils. Prenez un profilé en I, h = 300, b = 150, semelle t_f = 12, âme t_w = 8 mm, et décrivez-le de la seule façon qu'un programme puisse : comme les 12 points de coin de son contour, en traçant les semelles et les extrémités d'âme dans l'ordre. Pas de formule, pas de décomposition, juste le polygone.

Un passage de la famille du lacet sur ces 12 sommets renvoie :

• A = 5 808 mm²
• I_x = 88 709 184 mm⁴ ≈ 88,71 × 10⁶
• I_y = 6 761 776 mm⁴ ≈ 6,76 × 10⁶
• S_x = 591 395 mm³, r_x = 123,59 mm, r_y = 34,12 mm

Vérifions maintenant à la manière du manuel, le composite de trois rectangles avec le théorème des axes parallèles, I_x = [b·h³ − (b−t_w)·h_w³] / 12 avec h_w = 276 : [150·300³ − 142·276³] / 12 = 88 709 184 mm⁴. Les deux concordent à chaque chiffre, non pas à trois décimales, à chaque chiffre, parce que l'intégrale polygonale et l'intégrale composite sont la même intégrale évaluée de deux façons.

Cette exactitude est la raison pour laquelle une seule fonction polygonale peut piloter une base de profils entière. Elle se compose aussi : un PRS, une poutre alvéolaire à trous d'âme hexagonaux (voir poutres ajourées et alvéolaires), un caisson à vide interne, tous ne sont que plus de sommets, ou une boucle extérieure plus une boucle de trou. Le profilé en I à 12 sommets n'est pas un cas particulier dans le code ; c'est le cas général à 12 coins.

Tout ce qui précède était exact parce que les arêtes étaient droites. Un cercle, un congé, une paroi de tube a un contour courbe, et c'est le seul endroit où une vraie intégration numérique, au sens approché, entre réellement. Le geste universel est de remplacer la courbe par un polygone à nombreux côtés et d'exécuter dessus la formule polygonale exacte. La seule question est le nombre de côtés nécessaires, et la réponse est un résultat de convergence net.

Approchez un cercle de rayon 100 mm (A exact = πr² = 31 415,9 mm², I exact = πr⁴/4 = 78,54 × 10⁶ mm⁴) par un polygone régulier inscrit à n côtés :

• n = 6 : A écarté de −17,3 %, I écarté de −31,1 % (un hexagone est un mauvais cercle)
• n = 12 : A écarté de −4,51 %, I écarté de −8,77 %
• n = 24 : A écarté de −1,14 %, I écarté de −2,26 %
• n = 48 : A écarté de −0,285 %, I écarté de −0,570 %
• n = 96 : A écarté de −0,071 %, I écarté de −0,143 %
• n = 360 : A écarté de −0,005 %, I écarté de −0,010 %

Le motif est exact et mérite d'être assimilé : l'erreur décroît en 1/n². Passer de 12 à 24 côtés a réduit l'erreur d'aire de 4,51 % à 1,14 %, un facteur 3,96, en pratique le 4 théorique. Chaque doublement des côtés divise l'erreur par quatre. C'est une convergence du second ordre, le même ordre que celui rencontré dans la famille du trapèze et de Simpson en analyse, et c'est pourquoi un moteur de section peut atteindre cinq chiffres significatifs de précision sur un rond avec quelques centaines de côtés et n'y plus jamais penser.

L'alternative, et ce qu'un solveur de contraintes utilise sur l'intérieur de l'élément plutôt que sur le contour de section, est la quadrature de Gauss : au lieu de nombreux côtés égaux, échantillonner l'intégrande en quelques points placés de façon optimale. Comme y² n'est que quadratique, une règle de Gauss à 2 points l'intègre exactement sur un élément à côtés droits. Même idée d'analyse, l'échantillonnage pour remplacer une intégrale, réglée pour un autre travail.

Deux extensions couvrent presque toute section réelle, et les deux tombent tout droit de l'aire signée.

Trous par soustraction

Une section creuse, un tube, un caisson, une plaque à trou de boulon, est un contour extérieur moins un contour intérieur. Comme l'aire du lacet est signée, vous listez la boucle extérieure dans le sens trigonométrique (positive) et la boucle du trou dans le sens horaire (négative), vous additionnez chaque caractéristique des deux, et le trou se soustrait tout seul. Il n'y a pas de code spécial pour les trous : A = A_ext − A_trou, et la même annulation se produit pour les moments statiques et quadratiques avant l'unique transport au centre de gravité à la fin. Un tube rectangulaire est une boucle extérieure à quatre sommets et une boucle intérieure à quatre sommets, huit nombres, et le I exact d'un RHS en sort.

Composites par sommation

L'autre flux courant est l'élément reconstitué soudé : une poutre est une plaque de semelle supérieure, une plaque de semelle inférieure et une plaque d'âme. Ici vous n'avez même pas besoin de la voie polygonale ; le théorème des axes parallèles sur des rectangles suffit. Pour chaque plaque, ajoutez son propre b·h³/12 plus A·d² au total courant, où d est la distance du centre de gravité de la plaque à celui de la section. Sommez les aires, sommez les moments statiques pour trouver le centre de gravité, puis sommez les inerties transportées. C'est le même terme de Steiner qu'avant, désormais utilisé dans le sens direct, un rectangle à la fois.

La seule chose à surveiller : des plaques qui se chevauchent comptent deux fois. Si deux rectangles d'un composite partagent une aire intérieure, la sommation ajoute cette région deux fois et gonfle silencieusement A et I. Des arêtes qui se touchent, le cas soudé normal, sont acceptables ; un vrai chevauchement est une erreur de modélisation que le moteur signale plutôt que d'avaler.

Les mathématiques sont courtes, donc les erreurs sont toujours dans la comptabilité. Celles qui mordent, par ordre de fréquence d'apparition :

• Sens de parcours. Listez le polygone extérieur dans le sens horaire et l'aire signée sort négative. Un moteur robuste force la boucle extérieure positive et le trou négatif quelle que soit votre saisie, mais une implémentation à la main qui fait confiance au signe brut rapportera une aire négative et, pire, une inertie plausible mais fausse.
• Oublier le transport au centre de gravité. Les sommes en un passage donnent l'inertie par rapport à l'origine. Si vous vous arrêtez là et rapportez I_x0 comme I_x, vous avez surestimé la rigidité de A·C_y², ce qui pour une section décalée de l'origine peut dépasser la vraie réponse. La soustraction n'est pas optionnelle.
• Un polygone ouvert. La boucle suppose que le dernier sommet se reconnecte au premier. Fournissez-lui une liste non fermée, ou répétez par accident le sommet de fermeture sans que le code s'y attende, et une arête fantôme corrompt chaque somme. Les bons moteurs dédoublonnent un point de fermeture explicitement répété et ferment le reste implicitement.
• Contours auto-sécants. Un huit ou un contour qui se croise viole l'hypothèse de polygone simple ; les triangles signés ne partitionnent plus la région et le résultat n'a aucun sens. Tracez le contour une fois, sans croisement.
• Unités. L'inertie varie comme la longueur à la puissance quatre. Saisissez le contour en centimètres quand le reste du modèle est en millimètres et I est faux d'un facteur 10⁴. Gardez une seule unité de longueur d'un bout à l'autre ; un moteur de section qui travaille en interne en mm⁴ attend des mm en entrée.

Aucune de ces erreurs n'est subtile une fois nommée, et chacune a produit une section qui semblait correcte et se vérifiait fausse. L'intérêt d'utiliser le moteur livré est que les cinq sont gérées pour vous.

Tout ce qui précède est en service dans le calculateur de moment d'inertie de CalcSteel. Le constructeur de section personnalisée prend soit une liste de plaques soudées (la voie composite), soit les sommets bruts de n'importe quel polygone, éventuellement avec un trou, et renvoie A, le centre de gravité, I_x, I_y, I_xy, les axes principaux, S et r, dessinés à l'échelle, et confrontés en direct à la vraie base de 1 309 profils. C'est exactement la boucle du théorème de Green de cet article, la même fonction qui calcule les caractéristiques derrière chaque page de profil, si bien que le nombre que vous obtenez dans le navigateur est le nombre qu'utilise le moteur d'analyse.

Saisissez le rectangle 200 × 300 et regardez revenir 60 000 mm² et 450 × 10⁶ mm⁴. Collez les 12 sommets du profilé en I et lisez 88,71 × 10⁶. Tracez la cornière à ailes égales et voyez les axes principaux basculer à 45°. Il n'y a pas de maillage, pas d'attente, parce qu'il n'y a aucune approximation à faire converger, juste un passage sur les arêtes.

Une fois la section obtenue, l'étape suivante est ce que la structure en fait : injectez I dans une vérification de flèche de poutre ou la charge de flambement d'Euler, ou lisez comment le même I pilote les limites de flèche à l'état de service. Les caractéristiques de section sont le pont entre l'analyse et un élément qui passe la vérification.