Régime stationnaire

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 tutoriel/thermique.py ; tous les extraits ci-dessous en sont issus directement, dans l’ordre du fichier.

L’équation résolue

La page Thermique linéaire a posé l’équation de la chaleur et sa forme discrétisée, \( [C]\{\dot{T}\} + [K]\{T\} = \{P\} \). En régime stationnaire la température ne varie plus : le terme de capacité disparaît, et il reste l’équilibre entre la conduction et les sollicitations,

\[ \operatorname{div}(-\lambda \operatorname{grad} T) = q \ \text{ dans } V, \qquad T = T_{\text{imp}} \ \text{ sur } \partial V_T, \qquad -\phi \cdot n = \phi_{\text{imp}} + h\,(T_f - T) \ \text{ sur } \partial V_\phi \]

— la troisième équation compte la chaleur reçue par le bord (\( n \) la normale sortante) : un \( \phi_{\text{imp}} \) positif en apporte, les −40 kW/m² de la face gauche en retirent. 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 \]

La température imposée n’entre pas dans \( [K] \) : elle passe par des multiplicateurs de Lagrange \( \{\mu\} \), un par nœud bloqué. La relation \( [A]\{T\} = \{T_{\text{imp}}\} \) — chaque ligne de \( [A] \) lit la température d’un nœud de l’alésage — s’ajoute au système, et c’est celui-ci que pyrucast résout :

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

Les \( \{\mu\} \) sont les multiplicateurs, et leur opposé est la réaction : \( -[A]^T\{\mu\} \), les flux qu’il faut injecter aux nœuds bloqués pour les tenir à 250 °C.

Le second membre se vérifie par sa somme. Les fonctions de forme valent 1 à elles toutes en chaque point, si bien que la somme des lignes de \( \{P\} \) est la puissance que les sollicitations apportent à la pièce :

\[ \sum_i P_i = \int_V q \, dV + \int_{\partial V_\phi} (\phi_{\text{imp}} + h\,T_f) \, dS \]

— le contrôle que l’on fait plus bas sur le second membre de l’étape 2.

Chaque terme correspond à un objet du script :

Les trois dernières lignes disent deux choses. D’abord, un chargement est un terme du modèle : on le réunit au modèle, sa densité entre dans le champ matériau, et pc.node_field.external_forces calcule le second membre de tous les termes à la fois. Ensuite, flux est l’unique opérateur de charge répartie (première et dernière de ces lignes), quelle que soit la dimension du sous-espace éléments finis sur lequel on l’applique — une face QUA8 intègre une densité surfacique, un HEX20 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 tutoriel/maillage.py et calcule sur le volume HEX20 structuré de la page précédente — 640 hexaèdres.

    # FR — Le maillage de tutoriel/maillage.py, tel quel ; `plot=False` : pas ses figures.
    # EN — tutoriel/maillage.py'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 de la page Maillage.

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, 320 nœuds : 120 sommets et 200 nœuds milieux, puisque le volume est quadratique.

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 320 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 pc.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 QUA8 portés.
    # EN — The left face: the nodes of the plane x = 0, then the QUA8 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 QUA8 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é ; la contrainte crée ses multiplicateurs.
    # EN — Dirichlet: the constrained support; the constraint mints its multipliers.
    blocage_T = pc.model.dirichlet(alesage, modele, "T")
    modele = modele | blocage_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)]
    )

pc.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, c’est la relation \( [A]\{T\} = \{T_{\text{imp}}\} \) du système par blocs posé plus haut. pc.model.dirichlet ne reçoit qu’un maillage POI1 : le support bloqué (l’alésage, tel que points_on_cylinder l’a renvoyé). Les inconnues \( \{\mu\} \), elles, ont besoin de leurs propres nœuds : la contrainte les crée elle-même, un par nœud bloqué et à la même position, et blocage_T.multiplier_mesh() rend le maillage qui les porte. Deux jeux de nœuds distincts, donc deux jeux d’inconnues : la solution renvoyée contient les deux, la température « T » et les multiplicateurs « lambda », dont l’opposé est la réaction.

Le flux imposé, lui, est un terme du modèle : pc.model.flux le déclare sur l’espace éléments finis de la face gauche, et | le réunit au reste.

Le champ matériau se construit à partir du modèle : material_field sait quels coefficients celui-ci réclame — « k » pour la conduction, et la densité « phi_q » pour le flux.

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(blocage_T.multiplier_mesh(), ["imposed"])
    temperature_imposee[0].add_to_component("imposed", T_IMPOSEE)

Le flux imposé ne touche que le second membre : pc.node_field.external_forces intègre \( \int [N]^T \phi_{\text{imp}} \, dS \) sur les 20 QUA8 de la face gauche, avec la densité « phi_q » lue dans le champ matériau, et rend un champ nodal.

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 par blocs, en face des inconnues \( \{\mu\} \).

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 QUA8 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 HEX20 au lieu de QUA8. 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"])

    # FR — Le blocage de l'étape 1, repris tel quel avec ses multiplicateurs.
    # EN — Step 1's constraint, taken back as is with its multipliers.
    modele = modele | blocage_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 pc.model.boundary_transfer(basse_fes, conduction, ["T"]), 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, blocage_T, est repris tel quel, avec ses multiplicateurs : la température imposée posée sur eux à l’étape 1 reste donc valable.

Un seul material_field couvre le tout — « k » est réclamé par la conduction, « h_T » (le coefficient d’échange) et « a_ext_T » (la température du fluide) 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 HEX20, 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 HEX20 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 : le terme boundary_transfer le porte déjà, et le même appel pc.node_field.external_forces que pour le flux imposé le calcule, avec « h_T » et « a_ext_T » lus dans le champ matériau.

La source volumique, elle, est le terme \( \int_V [N]^T q \, dV \) : encore flux, mais appliqué à des HEX20. La dimension de l’intégrale est celle des éléments qu’on lui donne, donc une densité volumique ici. Un chargement ne contribue à aucune matrice : il forme ici un modèle à lui seul, source, avec son propre champ matériau, parce que sa densité « phi_q » n’est pas celle du flux de la face gauche.

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 ne sont jamais additionnés : chaque zone passe telle quelle, et quand deux zones donneraient deux valeurs différentes à un nœud partagé, l’opération s’arrête sur une erreur qui renvoie vers restrict_like. 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 Champ détaille cette algèbre ; la page Dilatation, matériau homogène 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, 64,8 W, vaut la somme des trois puissances prises séparément (−80 W de flux sortant, −115,2 W de convection, 260 W de source). L’union, elle, n’irait pas jusque-là : aux nœuds partagés les valeurs diffèrent, et elle s’arrêterait sur une erreur.

La température imposée, elle, vit sur le maillage des multiplicateurs, dont les nœuds, créés par la contrainte, n’appartiennent à aucun autre maillage : 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 de la page Maillage c’est la convection qui l’emporte : rien ne passe au-dessus de la température de l’alésage.

Ce régime établi est celui que la pièce atteint quand les sollicitations durent : la page Régime transitoire part d’une pièce à 0 °C, installe les mêmes chargements au fil du temps, et s’approche de ce champ au dernier instant.

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 :

"""Tutoriel débutant — 2. Régime stationnaire. / Beginner tutorial — 2. Steady state.

FR — Reprend la **chape percée** de la page « Maillage » — le volume HEX20
structuré est importé de `tutoriel/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`). Ce
script ne traite ni rayonnement ni terme transitoire : ce sont les scripts
`tutoriel/thermique_transitoire.py` et `tutoriel/thermique_non_lineaire.py`.
Le détail pas à pas est dans le livre, page « Régime stationnaire ».

EN — Picks the meshing page's **pierced lug** back up — the structured HEX20
volume is imported from `tutoriel/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`). This script handles
neither radiation nor a transient term: those are
`tutoriel/thermique_transitoire.py` and `tutoriel/thermique_non_lineaire.py`.
The step-by-step walkthrough lives in the book's steady-state page.

Lancement / Running ::

    pip install pyrucast
    python tutoriel/thermique.py

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

import os

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

# FR — La vue des figures de tutoriel/maillage.py : azimut, élévation, échelle.
# EN — The view of tutoriel/maillage.py'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 de tutoriel/maillage.py, tel quel ; `plot=False` : pas ses figures.
    # EN — tutoriel/maillage.py'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 QUA8 portés.
    # EN — The left face: the nodes of the plane x = 0, then the QUA8 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é ; la contrainte crée ses multiplicateurs.
    # EN — Dirichlet: the constrained support; the constraint mints its multipliers.
    blocage_T = pc.model.dirichlet(alesage, modele, "T")
    modele = modele | blocage_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(blocage_T.multiplier_mesh(), ["imposed"])
    temperature_imposee[0].add_to_component("imposed", 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"])

    # FR — Le blocage de l'étape 1, repris tel quel avec ses multiplicateurs.
    # EN — Step 1's constraint, taken back as is with its multipliers.
    modele = modele | blocage_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 HEX20, 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 HEX20 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} HEX20")
    print(f"face gauche    : {face_gauche.cell_count} QUA8")
    print(f"face convectée : {face_basse.cell_count} QUA8")
    print(f"zone chauffée  : {zone_source.cell_count} HEX20")
    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 : Régime transitoire.