Calcul thermique

Reprend la chape percée de la page Maillage pour un calcul de conduction thermique stationnaire, mené en deux temps : d’abord la conduction seule (température imposée, flux imposé), puis le même problème enrichi d’un film convectif et d’une source volumique. Chaque étape est tracée avant d’être résolue — les régions chargées d’abord, le champ de température ensuite.

Le script complet est formation/thermique.py ; tous les extraits ci-dessous en sont issus directement, dans l’ordre du fichier.

L’équation résolue

Le bilan d’énergie sur un volume V s’écrit, avec φ la densité de flux de chaleur et q une source volumique :

\[ \rho c_p \frac{\partial T}{\partial t} + \operatorname{div}(\phi) = q, \qquad \phi = -\lambda \operatorname{grad}(T) \]

En régime stationnaire le premier terme disparaît : il reste \( \operatorname{div}(-\lambda \operatorname{grad} T) = q \), complété par les conditions aux limites — température imposée sur une partie du bord, flux sur l’autre :

\[ T = T_{\text{imp}} \ \text{ sur } \partial V_T, \qquad \phi \cdot n = \phi_{\text{imp}} + h\,(T - T_f) \ \text{ sur } \partial V_\phi \]

Le terme en \( h \) est la convection (loi de Newton) : un échange avec un fluide à \( T_f \), proportionnel à l’écart de température. Il dépend de l’inconnue, donc il ne se range pas entièrement au second membre — on y revient plus bas.

Discrétisée sur les éléments finis, avec \( T(x) = [N(x)]\{T\} \) et \( \operatorname{grad}(T) = [B(x)]\{T\} \), l’équation devient un système linéaire :

\[ [K]\{T\} = \{P\} \]

\[ [K] = \int_V [B]^T \lambda [B] \, dV + \int_{\partial V_\phi} h\,[N]^T [N] \, dS \]

\[ \{P\} = \int_V [N]^T q \, dV + \int_{\partial V_\phi} [N]^T (\phi_{\text{imp}} + h\,T_f) \, dS \]

Chaque terme correspond à un objet du script :

Les trois dernières lignes disent le point notable : en pyrucast, flux est l’unique opérateur de charge répartie, quelle que soit la dimension du sous-espace éléments finis sur lequel on l’applique — une face QUA4 intègre une densité surfacique, un HEX8 une densité volumique. Flux surfacique, pression et source volumique passent donc tous par le même opérateur.

Les données du calcul

K_COND = 50.0  # W/m/K
FLUX_IMPOSE = -40_000.0  # W/m², face gauche / left face
H_CONV, T_EXT = 240.0, -80.0  # W/m²/K, °C — convection, face z = 0
SOURCE_VOLUMIQUE = 2600e3  # W/m³, zone chauffée / heated zone (≈ 260 W)
T_IMPOSEE = 250.0  # °C, alésage / bore

# FR — La zone chauffée est une tranche de la pièce, entre deux abscisses.
# EN — The heated zone is a slice of the part, between two abscissae.
SOURCE_X_MIN, SOURCE_X_MAX = 0.33 * LENGTH, 0.51 * LENGTH

# FR — Une face plane vaut zéro à l'arrondi près : on sélectionne une bande.
# EN — A flat face is zero up to rounding: a band is selected, not a value.
TOL = 1e-9  # m

Toutes les valeurs physiques sont groupées en tête de fichier, en unités SI : un acier (\( \lambda = 50 \) W/m/K), un flux sortant de −40 kW/m² sur la face gauche, un film convectif assez vif (\( h = 240 \) W/m²/K vers un fluide à −80 °C), une source de 2,6 MW/m³ dans la tranche chauffée (environ 260 W au total) et l’alésage tenu à 250 °C.

TOL mérite un mot : les nœuds d’une face plane valent zéro à l’arrondi machine près, jamais exactement zéro. Une sélection par coordonnée se fait donc toujours sur une bande, [-TOL, TOL], et la tolérance est ici explicite plutôt que cachée dans le mailleur.

On ne remaille pas : on importe

Le script ne redonne aucune cote. Il importe structured_mesh de formation/maillage.py et calcule sur le volume HEX8 structuré du chapitre précédent — 640 hexaèdres.

    # FR — Le maillage du chapitre 1, tel quel ; `plot=False` : pas ses figures.
    # EN — Chapter 1's mesh, as is; `plot=False`: without its figures.
    _, volume = structured_mesh(plot=False)

    # FR — Les charges réparties s'intègrent sur des faces : il faut la peau.
    # EN — Distributed loads integrate over faces: the skin is needed.
    peau = pc.mesh.consolidate(pc.mesh.skin(volume))

C’est la bonne façon d’enchaîner deux calculs sur une même pièce : un maillage reconstruit à l’identique dans deux scripts donnerait deux jeux de nœuds distincts, et toute condition posée sur l’un serait sans effet sur l’autre. plot=False demande seulement de ne pas retracer les figures du chapitre 1.

La deuxième ligne prépare la suite. Les chargements répartis s’intègrent sur des faces et non sur des nœuds — c’est le \( [N]^T \) des intégrales ci-dessus : il leur faut de vraies mailles de bord, que pyrucast.mesh.skin extrait du volume d’hexaèdres. pyrucast.mesh.consolidate ramène cette peau à un seul sous-maillage, pour que les sélections qui suivent en renvoient un seul elles aussi.

Étape 1 — conduction seule

Le premier problème n’a que deux conditions aux limites : l’alésage tenu à 250 °C, et un flux sortant de −40 kW/m² sur la face gauche. Aucun numéro de nœud n’apparaît dans le script : les régions sont découpées géométriquement, ici par forme, avec la famille pyrucast.mesh.points_*.

L’alésage, sur un cylindre

    # FR — L'axe du trou : la normale du plan de la pièce (Y), par le centre.
    # EN — The hole's axis: the part plane's normal (Y), through the centre.
    bas_axe = [LENGTH, -THICKNESS, HEIGHT / 2.0]
    haut_axe = [LENGTH, 2.0 * THICKNESS, HEIGHT / 2.0]

    # FR — L'alésage : les nœuds sur le cylindre, lus à même le volume.
    # EN — The bore: the nodes on the cylinder, read straight off the volume.
    alesage = pc.mesh.consolidate(
        pc.mesh.points_on_cylinder(volume, bas_axe, haut_axe, HOLE_RADIUS)
    )

L’axe du trou se donne par deux points, débordant de part et d’autre de la pièce : la normale du plan de la pièce (\( Y \)), passant par le centre du demi-disque. points_on_cylinder retient alors les nœuds sur le cylindre de rayon HOLE_RADIUS — les disques d’extrémité sont laissés de côté, ce sont des faces planes et points_on_plane est là pour celles-là. Ce qui revient est donc exactement la paroi du trou, 120 nœuds.

Un blocage ne demande que des nœuds. La sélection se lit donc directement sur le volume, sans en extraire la peau : les nœuds de la paroi du trou sont des nœuds de bord par définition, et points_on_cylinder y trouve les mêmes 120 nœuds que sur la peau.

Un points_* renvoie déjà un maillage POI1. Le résultat est utilisable tel quel comme support d’un model.dirichlet, sans passer par pyrucast.mesh.to_poi1. Seul mesh.consolidate reste nécessaire, pour écarter le sous-maillage vide que laisse la partie du volume qui ne touche pas le trou (le volume en compte deux : la grille et la couronne).

La face gauche, sur un plan

    # FR — La face gauche : les nœuds du plan x = 0, puis les QUA4 portés.
    # EN — The left face: the nodes of the plane x = 0, then the QUA4 they carry.
    noeuds_gauche = pc.mesh.points_on_plane(peau, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0])
    face_gauche = pc.mesh.elements_on(peau, noeuds_gauche, strict=True)

Même principe, mais un flux s’intègre sur une surface : les nœuds ne suffisent pas. pyrucast.mesh.elements_on(..., strict=True) remonte des nœuds sélectionnés aux mailles dont tous les sommets sont retenus — ici les 20 QUA4 de la face gauche. Le plan donné à points_on_plane est infini, mais il ne coupe la peau qu’à cet endroit.

La figure des conditions aux limites

    # FR — Une couleur par région, la peau en fil de fer autour.
    # EN — One colour per region, the skin drawn as a wireframe around them.
    alesage.unit().face_color = BLEU
    face_gauche.unit().face_color = ROUGE
    show(
        peau | alesage | face_gauche,
        "Étape 1 — conditions aux limites",
        "thermique-cl-conduction.svg",
        wireframe=True,
    )

Une couleur par région, et la figure devient le schéma. face_color se pose sur le sous-maillage (unit() le désigne quand il n’y en a qu’un), la peau est tracée en fil de fer autour (wireframe=True) : la figure ci-dessus se lit comme le croquis des conditions aux limites, sans annotation manuelle.

Le modèle et son matériau

    # FR — Le modèle porte « T » (primal) et « q » (dual) sur tout le volume.
    # EN — The model carries "T" (primal) and "q" (dual) over the whole volume.
    fes = pc.FiniteElementSpace(volume)
    modele = pc.model.heat_conduction(fes)

    # FR — Dirichlet : le support bloqué, et un jumeau pour les multiplicateurs.
    # EN — Dirichlet: the constrained support, and a twin for the multipliers.
    multiplicateur_T = pc.mesh.translate(alesage, [0.0, 0.0, 0.0])
    modele = modele | pc.model.dirichlet(modele, "T", alesage, multiplicateur_T)

    # FR — Le flux imposé sur la face gauche est un terme du modèle.
    # EN — The imposed flux on the left face is a term of the model.
    gauche_fes = pc.FiniteElementSpace(face_gauche)
    modele = modele | pc.model.flux(gauche_fes, modele, "q")

    # FR — La conduction réclame « k », la charge sa densité.
    # EN — Conduction asks for "k", the load for its density.
    materiaux = pc.element_field.material_field(
        modele, [("k", K_COND), ("phi_q", FLUX_IMPOSE)]
    )

model.heat_conduction déclare le couple de degrés de liberté « T » (primal) et « q » (dual) sur tout le volume et porte le terme \( \int_V [B]^T \lambda [B] \, dV \) ; pc.matrix.stiffness l’intègre réellement, avec le \( \lambda \) lu dans le champ matériau sous le nom « k ».

La température imposée passe par des multiplicateurs de Lagrange : le système résolu n’est plus \( [K]\{T\} = \{P\} \) mais

\[ \begin{bmatrix} K & C^T \\ C & 0 \end{bmatrix} \begin{Bmatrix} T \\ \lambda \end{Bmatrix} = \begin{Bmatrix} P \\ T_{\text{imp}} \end{Bmatrix} \]

où \( C \) est la relation \( T = T_{\text{imp}} \) sur les nœuds de l’alésage. D’où les deux maillages POI1 donnés à model.dirichlet : le support bloqué (l’alésage, tel que points_on_cylinder l’a renvoyé) et un jumeau qui porte les inconnues \( \lambda \), obtenu par copie translatée de zéro — deux jeux de nœuds distincts, donc deux jeux d’inconnues. La solution renvoyée contient les deux, et les multiplicateurs sont les réactions (ici les flux qu’il faut injecter pour tenir l’alésage à 250 °C).

Le champ matériau, lui, se construit à partir du modèle : material_field sait quels coefficients celui-ci réclame, et la conduction n’en demande qu’un, « k ».

Les deux chargements

    flux_gauche = pc.node_field.external_forces(modele, materiaux)

    # FR — Température imposée, posée sur le maillage des multiplicateurs.
    # EN — Imposed temperature, set on the multipliers' mesh.
    temperature_imposee = pc.NodeField(multiplicateur_T, ["imposed_T"])
    temperature_imposee[0].add_to_component("imposed_T", T_IMPOSEE)

Le flux imposé est un pur second membre : pc.model.flux intègre \( \int [N]^T \phi_{\text{imp}} \, dS \) sur les 20 QUA4 de la face gauche et rend un champ nodal. L’espace éléments finis se construit sur le sous-maillage de la face, et gauche_fes[0] en désigne l’unique zone.

La température imposée, elle, se pose sur le maillage des multiplicateurs et non sur l’alésage lui-même : c’est le \( T_{\text{imp}} \) du second bloc du système ci-dessus, en face des inconnues \( \lambda \).

Résolution

    # FR — `[K]{T} = {P}` : matrice assemblée, second membre réuni par `|`.
    # EN — `[K]{T} = {P}`: assembled matrix, right-hand side gathered by `|`.
    K = pc.matrix.stiffness(modele, materiaux)
    t_conduction = pc.solver.solve(K, flux_gauche | temperature_imposee)
    show_nodefield(
        volume, t_conduction, "Étape 1 — température (°C)", "thermique-conduction.svg"
    )

pc.matrix.stiffness intègre la matrice, | réunit les deux chargements — ils vivent sur des maillages disjoints, il n’y a donc rien à sommer — et pyrucast.solver.solve factorise la matrice creuse (LU parallèle) en mettant la factorisation en cache : deux résolutions sur la même matrice ne la factorisent qu’une fois.

Le résultat est le gradient attendu : 250 °C tenus à l’alésage, 32 °C sur la face gauche d’où la chaleur s’échappe.

Étape 2 — convection et source volumique

On ajoute maintenant les deux sollicitations restantes, sans rien retoucher aux précédentes : un film convectif sous la pièce, sur la face \( z = 0 \), et une tranche chauffée entre \( 0{,}33\,L \) et \( 0{,}51\,L \). Ces deux régions n’ont pas de forme simple à nommer : elles sont repérées par coordonnée, la seconde façon de découper une région.

La surface convectée, par coordonnée

    # FR — La face convectée, z = 0 : repérée par coordonnée, pas par forme.
    # EN — The convected face, z = 0: located by coordinate, not by shape.
    z_peau = pc.node_field.positions(peau, ["Z"])
    noeuds_bas = pc.mesh.select(z_peau, ge=-TOL, le=TOL)
    face_basse = pc.mesh.elements_on(peau, noeuds_bas, strict=True)
    face_basse.unit().face_color = TURQUOISE
    show(
        peau | face_basse,
        "Étape 2 — surface convectée",
        "thermique-cl-convection.svg",
        wireframe=True,
    )

Une coordonnée est un champ nodal comme un autre. pyrucast.node_field.positions(peau, ["Z"]) rend la coordonnée Z des nœuds de la peau sous forme de NodeField, et pyrucast.mesh.select garde ceux dont la valeur tombe dans une bande — ge=-TOL, le=TOL pour « z = 0 ». Le résultat est un maillage POI1, exactement comme celui d’un points_* : la suite ne change pas, elements_on(..., strict=True) remonte aux QUA4 que ces nœuds portent entièrement.

La zone chauffée, en bande

    # FR — La zone chauffée : même démarche sur X, en bande, et sur le volume.
    # EN — The heated zone: same approach on X, as a band, over the volume.
    x_volume = pc.node_field.positions(volume, ["X"])
    noeuds_source = pc.mesh.select(x_volume, ge=SOURCE_X_MIN, le=SOURCE_X_MAX)
    zone_source = pc.mesh.consolidate(
        pc.mesh.elements_on(volume, noeuds_source, strict=True)
    )
    zone_source.unit().face_color = VERT
    show(
        peau | zone_source,
        "Étape 2 — zone chauffée",
        "thermique-cl-source.svg",
        wireframe=True,
    )

Même démarche, mais sur X et sur le volume : une bande de valeurs au lieu d’une égalité, et des HEX8 au lieu de QUA4. Comme pour l’alésage, mesh.consolidate écarte le sous-maillage vide laissé par la partie du volume qui ne rencontre pas la bande.

strict=True approche la région par un escalier. La bande en X coupe le maillage entre deux abscisses quelconques, mais ce qui est retenu est le paquet des 80 hexaèdres dont tous les nœuds sont dedans : la tranche s’arrête donc aux frontières des éléments, bien visible sur la figure. C’est le prix à payer pour que la région chargée soit un sous-maillage conforme.

Le modèle complet

    # FR — La convection s'ajoute dans la matrice : `|` sur les mêmes ddl.
    # EN — Convection adds into the matrix: `|` on the very same dofs.
    basse_fes = pc.FiniteElementSpace(face_basse)
    conduction = pc.model.heat_conduction(fes)
    modele = conduction | pc.model.boundary_transfer(
        basse_fes, conduction, [("T", "q")]
    )
    modele = modele | pc.model.dirichlet(modele, "T", alesage, multiplicateur_T)

    # FR — Un seul champ matériau : « k » pour la conduction, « h » et son
    #      ambiant pour le film.
    # EN — A single material field: "k" for conduction, "h" and its ambient for
    #      the film.
    materiaux = pc.element_field.material_field(
        modele, [("k", K_COND), ("h_T", H_CONV), ("a_ext_T", T_EXT)]
    )

La convection est la seule des quatre sollicitations à toucher les deux membres du système, parce que \( \phi \cdot n = h\,(T - T_f) \) dépend de l’inconnue. Sa part en \( T \) donne \( \int h\,[N]^T[N] \, dS \), qui s’ajoute dans la matrice : ce n’est pas un système séparé, d’où le model.boundary_transfer(basse_fes, conduction, [("T", "q")]), construit contre la conduction puis réuni à elle par |, sur les mêmes degrés de liberté « T » et « q ». Le blocage de l’étape 1 est repris tel quel, avec les mêmes deux maillages POI1.

Un seul material_field couvre le tout — « k » est réclamé par la conduction, « h » par la convection.

Les deux nouveaux chargements

    # FR — Terme externe de la convection, h·T_ext : le modèle le porte.
    # EN — The convection's external term, h·T_ext: the model carries it.
    charge_convection = pc.node_field.external_forces(modele, materiaux)

    # FR — Source volumique sur des HEX8, donc une densité volumique. Une
    #      charge ne contribuant à aucune matrice, elle se tient très bien en
    #      modèle à elle seule — avec sa propre densité, distincte de celle du
    #      flux de bord bien qu'elles alimentent la même ligne « q ».
    # EN — A volume source over HEX8 cells. A load contributes to no matrix, so
    #      it stands perfectly well as a model of its own — with its own
    #      density, distinct from the boundary flux's though both feed "q".
    source_fes = pc.FiniteElementSpace(zone_source)
    source = pc.model.flux(source_fes, modele, "q")
    densite_source = pc.element_field.material_field(
        source, [("phi_q", SOURCE_VOLUMIQUE)]
    )
    charge_source = pc.node_field.external_forces(source, densite_source)

La part en \( T_f \) de la convection donne \( \int h\,T_f\,[N]^T \, dS \), un second membre ordinaire : c’est le même opérateur flux que pour le flux imposé, sur la surface convectée.

La source volumique, elle, est le terme \( \int_V [N]^T q \, dV \) : encore flux, mais appliqué à des HEX8. La dimension de l’intégrale est celle des éléments qu’on lui donne, donc une densité volumique ici.

Trois chargements qui se touchent : le second membre se somme

Le bas de la face gauche est sur \( z = 0 \), et la tranche chauffée débouche elle aussi sous la pièce : les trois chargements répartis partagent des nœuds. Leurs contributions doivent donc s’additionner là — et c’est précisément ce que l’union | ne fait pas.

Piège : deux régions chargées adjacentes. Leurs contributions nodales ne sont pas sommées automatiquement à l’assemblage : chaque chargement est assemblé sur son propre support, et l’union (|) juxtapose ces supports sans les additionner. À un nœud partagé, le solveur lit le second membre zone par zone et retient la valeur de la première qui définit le couple (nœud, composante) — l’autre contribution est perdue. L’union ne lève une erreur que si les deux valeurs diffèrent ; quand elles coïncident, elle passe sans rien dire.

Sommer les champs (+) ne suffit pas non plus tel quel : l’arithmétique de champs apparie elle aussi les zones par support (deux supports distincts sont recopiés tels quels), et elle ne fait même pas la vérification de cohérence de |. Il faut d’abord ramener les champs sur un support commun — pyrucast.node_field.restrict sur un même maillage retombe sur le support POI1 canonique de ce maillage, donc restrict(a, m) + restrict(b, m) est bien une somme nœud à nœud (la page Champs détaille cette algèbre ; le chapitre 3 en donne un exemple avec restrict_like).

    # FR — Les trois charges se touchent : support commun, puis `+` somme.
    # EN — The three loads touch: a common support first, then `+` really sums.
    noeuds_charges = pc.mesh.consolidate(
        pc.mesh.to_poi1(face_gauche | face_basse | zone_source)
    )
    second_membre = (
        pc.node_field.restrict(flux_gauche, noeuds_charges)
        + pc.node_field.restrict(charge_convection, noeuds_charges)
        + pc.node_field.restrict(charge_source, noeuds_charges)
    ) | temperature_imposee

D’où ces quelques lignes : le maillage POI1 de tous les nœuds chargés (to_poi1 puis mesh.consolidate, pour n’avoir qu’un seul sous-maillage donc un seul support), les trois champs restreints dessus, et leur somme par +. La vérification est immédiate — la puissance totale du second membre vaut la somme des trois puissances prises séparément, ce que l’union perdrait.

La température imposée, elle, vit sur le maillage des multiplicateurs, translaté donc disjoint de tous les autres : elle se réunit au reste par | sans rien avoir à sommer.

Résolution

    # FR — Même schéma qu'à l'étape 1, sur la matrice enrichie du terme convectif.
    # EN — Same pattern as step 1, on the matrix enriched with the film term.
    K = pc.matrix.stiffness(modele, materiaux)
    t_complet = pc.solver.solve(K, second_membre)
    show_nodefield(
        volume, t_complet, "Étape 2 — température (°C)", "thermique-complet.svg"
    )

Rien de nouveau par rapport à l’étape 1 : la même paire stiffness / solve, sur une matrice qui porte en plus le terme convectif.

La lecture du résultat suit les quatre chargements : l’alésage est toujours tenu à 250 °C par le blocage et la face gauche reste le point froid sous le flux sortant, mais les isothermes, rectilignes à l’étape 1, s’infléchissent maintenant autour de la tranche chauffée. Les deux nouvelles sollicitations jouent en sens contraire — la source réchauffe la moitié gauche (le minimum remonte de 32 à 40 °C), le film convectif pompe la chaleur sous la pièce — et avec les cotes du chapitre 1 c’est la convection qui l’emporte : rien ne passe au-dessus de la température de l’alésage.

Non disponible dans pyrucast.

• Rayonnement. Pas de condition de bord de type ϕ·n = εσ(T⁴∞ − T⁴) — seules conduction et convection (film/Robin) existent.
• Régime transitoire. La matrice de capacité \( [C] = \int_V \rho c_p [N]^T [N] \, dV \) est assemblable (pyrucast.matrix.mass), mais rien ne la relie encore à une boucle en temps : chaque pas résout \( [K]\{T\} = \{P\} \) stationnaire (pas de \( [C]\{\dot{T}\} + [K]\{T\} = \{P\} \) intégré en temps). Un Evolution peut faire varier un chargement stationnaire d’un pas à l’autre (voir Calcul mécanique pour ce mécanisme appliqué à la mécanique), mais c’est une suite de problèmes stationnaires indépendants, pas une intégration temporelle.

Script complet

Une fois les explications retirées, tout tient en une page — c’est le déroulé complet, du maillage importé aux deux résolutions :

"""Formation débutant — 2. Calcul thermique. / Beginner training — 2. Thermal.

FR — Reprend la **chape percée** du chapitre 1 — le volume HEX8 structuré est
importé de `formation/maillage.py`, aucune cote n'est redonnée — et y résout la
conduction **stationnaire** `div(-k·grad T) = q`, en deux temps :

1. **conduction seule** — température imposée sur l'alésage, flux imposé sur la
   face gauche ;
2. **conduction + convection + source** — film convectif sous la pièce et
   tranche chauffée, sans rien retoucher au reste.

Chaque étape est tracée avant d'être résolue, et les régions chargées sont
repérées **par leur géométrie** : par forme (`pyrucast.mesh.points_*`) ou par
coordonnée (`pyrucast.node_field.positions` + `pyrucast.mesh.select`). Ni
rayonnement ni terme transitoire. Le détail pas à pas est dans le livre, page
« Calcul thermique ».

EN — Picks chapter 1's **pierced lug** back up — the structured HEX8 volume is
imported from `formation/maillage.py`, not one dimension is restated — and
solves **steady** conduction `div(-k·grad T) = q` on it, in two steps:

1. **conduction alone** — imposed temperature on the bore, imposed flux on the
   left face;
2. **conduction + convection + source** — convective film under the part and
   heated slice, nothing else changes.

Each step is plotted before being solved, and the loaded regions are located
**by geometry**: by shape (`pyrucast.mesh.points_*`) or by coordinate
(`pyrucast.node_field.positions` + `pyrucast.mesh.select`). No radiation, no
transient term. The step-by-step walkthrough lives in the book's thermal page.

Lancement / Running ::

    maturin develop --release
    python formation/thermique.py

    # Figures du livre / book figures (book/src/formation/img/) :
    # PYRUCAST_FORMATION_IMG_DIR=book/src/formation/img python formation/thermique.py
"""

import os

import pyrucast as pc
from maillage import HEIGHT, HOLE_RADIUS, LENGTH, OUT, THICKNESS, show, structured_mesh

# FR — Vue commune aux figures du chapitre 1 : azimut, élévation, échelle.
# EN — The view shared by chapter 1's figures: azimuth, elevation, scale.
VUE = (-45, 25, 1.0)

# FR — Une couleur par région chargée, tenue d'une figure à l'autre.
# EN — One colour per loaded region, kept from one figure to the next.
BLEU = (0, 0, 255)  # alésage / bore
ROUGE = (255, 0, 0)  # face gauche / left face
TURQUOISE = (0, 190, 190)  # face convectée / convected face
VERT = (0, 170, 0)  # zone chauffée / heated zone

# ── Données physiques / Physical data ──────────────────────────────────────
K_COND = 50.0  # W/m/K
FLUX_IMPOSE = -40_000.0  # W/m², face gauche / left face
H_CONV, T_EXT = 240.0, -80.0  # W/m²/K, °C — convection, face z = 0
SOURCE_VOLUMIQUE = 2600e3  # W/m³, zone chauffée / heated zone (≈ 260 W)
T_IMPOSEE = 250.0  # °C, alésage / bore

# FR — La zone chauffée est une tranche de la pièce, entre deux abscisses.
# EN — The heated zone is a slice of the part, between two abscissae.
SOURCE_X_MIN, SOURCE_X_MAX = 0.33 * LENGTH, 0.51 * LENGTH

# FR — Une face plane vaut zéro à l'arrondi près : on sélectionne une bande.
# EN — A flat face is zero up to rounding: a band is selected, not a value.
TOL = 1e-9  # m


def show_nodefield(mesh: pc.Mesh, field: pc.NodeField, title: str, file: str) -> None:
    """FR — Trace `field` sur `mesh` : fenêtre interactive, ou SVG si `OUT`.

    EN — Plot `field` over `mesh`: interactive window, or SVG when `OUT` is set.
    """
    mesh.plot(
        view=VUE,
        title=title,
        field=field,
        component="T",
        cmap="viridis",
        smooth=1,
        save=os.path.join(OUT, file) if OUT else None,
    )


def main() -> None:
    # FR — Le maillage du chapitre 1, tel quel ; `plot=False` : pas ses figures.
    # EN — Chapter 1's mesh, as is; `plot=False`: without its figures.
    _, volume = structured_mesh(plot=False)

    # FR — Les charges réparties s'intègrent sur des faces : il faut la peau.
    # EN — Distributed loads integrate over faces: the skin is needed.
    peau = pc.mesh.consolidate(pc.mesh.skin(volume))

    # ── Étape 1 : régions / Step 1: regions ─────────────────────────────────
    # FR — L'axe du trou : la normale du plan de la pièce (Y), par le centre.
    # EN — The hole's axis: the part plane's normal (Y), through the centre.
    bas_axe = [LENGTH, -THICKNESS, HEIGHT / 2.0]
    haut_axe = [LENGTH, 2.0 * THICKNESS, HEIGHT / 2.0]

    # FR — L'alésage : les nœuds sur le cylindre, lus à même le volume.
    # EN — The bore: the nodes on the cylinder, read straight off the volume.
    alesage = pc.mesh.consolidate(
        pc.mesh.points_on_cylinder(volume, bas_axe, haut_axe, HOLE_RADIUS)
    )

    # FR — La face gauche : les nœuds du plan x = 0, puis les QUA4 portés.
    # EN — The left face: the nodes of the plane x = 0, then the QUA4 they carry.
    noeuds_gauche = pc.mesh.points_on_plane(peau, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0])
    face_gauche = pc.mesh.elements_on(peau, noeuds_gauche, strict=True)

    # FR — Une couleur par région, la peau en fil de fer autour.
    # EN — One colour per region, the skin drawn as a wireframe around them.
    alesage.unit().face_color = BLEU
    face_gauche.unit().face_color = ROUGE
    show(
        peau | alesage | face_gauche,
        "Étape 1 — conditions aux limites",
        "thermique-cl-conduction.svg",
        wireframe=True,
    )

    # ── Étape 1 : calcul / Step 1: analysis ─────────────────────────────────
    # FR — Le modèle porte « T » (primal) et « q » (dual) sur tout le volume.
    # EN — The model carries "T" (primal) and "q" (dual) over the whole volume.
    fes = pc.FiniteElementSpace(volume)
    modele = pc.model.heat_conduction(fes)

    # FR — Dirichlet : le support bloqué, et un jumeau pour les multiplicateurs.
    # EN — Dirichlet: the constrained support, and a twin for the multipliers.
    multiplicateur_T = pc.mesh.translate(alesage, [0.0, 0.0, 0.0])
    modele = modele | pc.model.dirichlet(modele, "T", alesage, multiplicateur_T)

    # FR — Le flux imposé sur la face gauche est un terme du modèle.
    # EN — The imposed flux on the left face is a term of the model.
    gauche_fes = pc.FiniteElementSpace(face_gauche)
    modele = modele | pc.model.flux(gauche_fes, modele, "q")

    # FR — La conduction réclame « k », la charge sa densité.
    # EN — Conduction asks for "k", the load for its density.
    materiaux = pc.element_field.material_field(
        modele, [("k", K_COND), ("phi_q", FLUX_IMPOSE)]
    )

    flux_gauche = pc.node_field.external_forces(modele, materiaux)

    # FR — Température imposée, posée sur le maillage des multiplicateurs.
    # EN — Imposed temperature, set on the multipliers' mesh.
    temperature_imposee = pc.NodeField(multiplicateur_T, ["imposed_T"])
    temperature_imposee[0].add_to_component("imposed_T", T_IMPOSEE)

    # FR — `[K]{T} = {P}` : matrice assemblée, second membre réuni par `|`.
    # EN — `[K]{T} = {P}`: assembled matrix, right-hand side gathered by `|`.
    K = pc.matrix.stiffness(modele, materiaux)
    t_conduction = pc.solver.solve(K, flux_gauche | temperature_imposee)
    show_nodefield(
        volume, t_conduction, "Étape 1 — température (°C)", "thermique-conduction.svg"
    )

    # ── Étape 2 : régions / Step 2: regions ─────────────────────────────────
    # FR — La face convectée, z = 0 : repérée par coordonnée, pas par forme.
    # EN — The convected face, z = 0: located by coordinate, not by shape.
    z_peau = pc.node_field.positions(peau, ["Z"])
    noeuds_bas = pc.mesh.select(z_peau, ge=-TOL, le=TOL)
    face_basse = pc.mesh.elements_on(peau, noeuds_bas, strict=True)
    face_basse.unit().face_color = TURQUOISE
    show(
        peau | face_basse,
        "Étape 2 — surface convectée",
        "thermique-cl-convection.svg",
        wireframe=True,
    )

    # FR — La zone chauffée : même démarche sur X, en bande, et sur le volume.
    # EN — The heated zone: same approach on X, as a band, over the volume.
    x_volume = pc.node_field.positions(volume, ["X"])
    noeuds_source = pc.mesh.select(x_volume, ge=SOURCE_X_MIN, le=SOURCE_X_MAX)
    zone_source = pc.mesh.consolidate(
        pc.mesh.elements_on(volume, noeuds_source, strict=True)
    )
    zone_source.unit().face_color = VERT
    show(
        peau | zone_source,
        "Étape 2 — zone chauffée",
        "thermique-cl-source.svg",
        wireframe=True,
    )

    # ── Étape 2 : calcul / Step 2: analysis ─────────────────────────────────
    # FR — La convection s'ajoute dans la matrice : `|` sur les mêmes ddl.
    # EN — Convection adds into the matrix: `|` on the very same dofs.
    basse_fes = pc.FiniteElementSpace(face_basse)
    conduction = pc.model.heat_conduction(fes)
    modele = conduction | pc.model.boundary_transfer(
        basse_fes, conduction, [("T", "q")]
    )
    modele = modele | pc.model.dirichlet(modele, "T", alesage, multiplicateur_T)

    # FR — Un seul champ matériau : « k » pour la conduction, « h » et son
    #      ambiant pour le film.
    # EN — A single material field: "k" for conduction, "h" and its ambient for
    #      the film.
    materiaux = pc.element_field.material_field(
        modele, [("k", K_COND), ("h_T", H_CONV), ("a_ext_T", T_EXT)]
    )

    # FR — Terme externe de la convection, h·T_ext : le modèle le porte.
    # EN — The convection's external term, h·T_ext: the model carries it.
    charge_convection = pc.node_field.external_forces(modele, materiaux)

    # FR — Source volumique sur des HEX8, donc une densité volumique. Une
    #      charge ne contribuant à aucune matrice, elle se tient très bien en
    #      modèle à elle seule — avec sa propre densité, distincte de celle du
    #      flux de bord bien qu'elles alimentent la même ligne « q ».
    # EN — A volume source over HEX8 cells. A load contributes to no matrix, so
    #      it stands perfectly well as a model of its own — with its own
    #      density, distinct from the boundary flux's though both feed "q".
    source_fes = pc.FiniteElementSpace(zone_source)
    source = pc.model.flux(source_fes, modele, "q")
    densite_source = pc.element_field.material_field(
        source, [("phi_q", SOURCE_VOLUMIQUE)]
    )
    charge_source = pc.node_field.external_forces(source, densite_source)

    # FR — Les trois charges se touchent : support commun, puis `+` somme.
    # EN — The three loads touch: a common support first, then `+` really sums.
    noeuds_charges = pc.mesh.consolidate(
        pc.mesh.to_poi1(face_gauche | face_basse | zone_source)
    )
    second_membre = (
        pc.node_field.restrict(flux_gauche, noeuds_charges)
        + pc.node_field.restrict(charge_convection, noeuds_charges)
        + pc.node_field.restrict(charge_source, noeuds_charges)
    ) | temperature_imposee

    # FR — Même schéma qu'à l'étape 1, sur la matrice enrichie du terme convectif.
    # EN — Same pattern as step 1, on the matrix enriched with the film term.
    K = pc.matrix.stiffness(modele, materiaux)
    t_complet = pc.solver.solve(K, second_membre)
    show_nodefield(
        volume, t_complet, "Étape 2 — température (°C)", "thermique-complet.svg"
    )

    print(f"volume         : {volume.cell_count()} HEX8")
    print(f"face gauche    : {face_gauche.cell_count()} QUA4")
    print(f"face convectée : {face_basse.cell_count()} QUA4")
    print(f"zone chauffée  : {zone_source.cell_count()} HEX8")
    print(
        f"étape 1 : T min = {t_conduction.min('T'):.1f} °C, "
        f"T max = {t_conduction.max('T'):.1f} °C"
    )
    print(
        f"étape 2 : T min = {t_complet.min('T'):.1f} °C, "
        f"T max = {t_complet.max('T'):.1f} °C"
    )


if __name__ == "__main__":
    main()

Suite : Calcul mécanique, qui réutilise ce champ de température.