Conduction thermique

Cette page décrit la physique de conduction thermique (HeatConduction) et la convection de surface associée. Elle suit le plan standard des physiques puis la déroule sur un exemple complet — une ligne chauffée par une source à une extrémité et maintenue à température fixe à l’autre — comparé à la solution analytique.

Pour la mécanique générique du Model (orchestration, DOFs, assemblage), voir Modèle physique. Ici on se concentre sur le cas thermique.

Équations continues résolues

En régime stationnaire, la forme forte est

\[ -\nabla\cdot\big(k\,\nabla T\big) = 0, \]

et en régime transitoire, l’équation de la chaleur porte un terme de stockage :

\[ \rho\,c_p\,\frac{\partial T}{\partial t} - \nabla\cdot(k\,\nabla T) = Q. \]

La forme variationnelle de Galerkine (multiplication par une température virtuelle, intégration par parties) fait apparaître la rigidité (conduction) et, en transitoire, la capacité (stockage) — leurs formes discrètes ci-dessous.

Forme discrétisée

La conductivité donne, cellule par cellule, la matrice de rigidité :

\[ K_{ij} = \int_K k(x)\,\nabla N_i\cdot\nabla N_j\,dx \quad\approx\quad \sum_g k(\xi_g)\,(\nabla N_i\cdot\nabla N_j)\big|_g\,|J|_g\,w_g \]

(implémentée dans src/models/heat_conduction.rs). En notant \( B = [\nabla N_1, \dots, \nabla N_n] \) la matrice des gradients de forme (taille \( d\times n \)), on a aussi \( K = \int_\Omega k\, B^\top B\, d\Omega \). Le bloc local est écrit aux positions row = (NodeId_i, "q"), col = (NodeId_j, "T"). Pour un SEG2 de longueur \(L\) et \(k\) uniforme on retrouve la matrice analytique \((k/L)\,[[1,-1],[-1,1]]\).

En transitoire, le terme de stockage discrétise en une matrice de capacité (l’analogue thermique de la matrice de masse, Cast3M CAPA) :

\[ C_{ij} = \int_\Omega \rho\,c_p\,N_i\,N_j\,d\Omega \;\approx\; \sum_g \rho\,c_p\,N_i(\xi_g)\,N_j(\xi_g)\,|J|_g\,w_g, \]

assemblée par assemble.mass (matériau rho, cp) et concentrable en diagonale par lump. Le système semi-discret est \( C\,\dot T + K\,T = F \) ; l’intégration en temps (θ-schéma, Euler implicite (C/\Delta t + K)) se pilote dans la couche Python.

Variables et matériau

La conductivité peut être orientée — voir Conduction orthotrope et anisotrope plus bas ; "k" est alors remplacée par les constantes de la symétrie choisie.

Mise en donnée (Rust, testé)

Le pipeline est toujours le même :

1. Coords — l’espace des nœuds (dimension géométrique).
2. Mesh — les éléments (ici des SEG2 alignés sur \([0,1]\)).
3. FiniteElementSpace — l’interpolation (lagrange1).
4. Matériau — un ElementField portant la composante "k", fabriqué commodément par element_field::material_field(&model, &[("k", …)]) (les sous-modèles sans matériau, comme Dirichlet, sont ignorés).
5. Model — model::heat_conduction(&fes), composé par | (union) avec les conditions limites.
6. Conditions limites :
   • Dirichlet (T imposée) : un sous-modèle model::dirichlet qui impose la valeur via multiplicateurs de Lagrange. L’utilisateur fournit le maillage des nœuds imposés et le maillage support des multiplicateurs — typiquement fabriqué depuis le premier avec le mesher générique barycenter (nœuds neufs colocalisés). La valeur imposée \(u_d\) s’écrit dans le chargement au slot imposed_T du nœud-multiplicateur (cf. Modèle physique).
   • Neumann / source : une charge ponctuelle est une valeur du chargement sur la composante duale "q" au nœud concerné ; un flux réparti sur un bord (ou un volume) se transforme en charges nodales cohérentes par l’opérateur flux (analogue de FLUX/PRES de Cast3M).
7. Assemblage + résolution — matrix::stiffness puis le solveur solver::lu::solve (LU creuse directe, voir Modèle physique).

Exemple : ligne chauffée

Problème. Sur \([0,1]\), une source de chaleur (flux de Neumann \(Q\)) est appliquée en \(x=0\), et la température est imposée à \(T=20\) en \(x=1\).

Solution analytique. Sans génération volumique, \(T’’=0\) : le profil est linéaire. En notant \(Q\) le flux injecté et \(k\) la conductivité,

\[ u(x) = 20 + \frac{Q}{k}\,(1 - x). \]

De plus, le multiplicateur de Lagrange au nœud imposé (la réaction qui maintient \(T=20\)) vaut exactement \(Q\) : tout le flux injecté en \(x=0\) ressort en \(x=1\) — un bilan d’énergie discret.

Code. L’exemple ci-dessous est le test d’intégration tests/thermal_line.rs : il est compilé et exécuté à chaque cargo test, donc garanti à jour avec l’API.

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn thermal_line_recovers_analytical_solution() -> Result<()> {
    // ── Problem data ───────────────────────────────────────────────────────
    const K: f64 = 1.0; // conductivité
    const Q: f64 = 10.0; // source de chaleur (flux de Neumann) en x = 0
    const T_IMPOSED: f64 = 20.0; // température imposée en x = 1
    const N_ELEMS: usize = 4;
    let h = 1.0 / N_ELEMS as f64;

    // ── Mesh: a line of SEG2 on [0, 1] ─────────────────────────────────────
    let coords = Handle::new(Coords::new(1)?);
    let nodes: Vec<Node> = (0..=N_ELEMS)
        .map(|i| Node::create_in(coords.clone(), &[i as f64 * h]))
        .collect::<Result<_>>()?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for i in 0..N_ELEMS {
        mesh.add_cell(&[nodes[i].id(), nodes[i + 1].id()])?;
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : conduction + Dirichlet T = 20 en x = 1 ────────────────────
    // The multipliers' support is built from the imposed node by the
    // `barycenter` mesher (a fresh co-located node). The model creates nothing.
    let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(
        nodes.last().unwrap(),
    ))?);
    let multiplier = mesh::barycenter(&imposed)?;
    let mult = multiplier.node(0, 0, 0)?.id();

    let conduction = model::heat_conduction(&fes)?;
    let dirichlet = model::dirichlet(&conduction, "T", &imposed, &multiplier, Default::default())?;
    let model = conduction.union(&dirichlet)?;

    // ── Material: uniform k (Dirichlet is skipped automatically) ───────────
    let materials = pyrucast::ops::element_field::material_field(&model, &[("k", K)])?;

    // ── Chargement : source Q en x = 0 (composante duale "q"), valeur imposée
    //    T = 20 at the multiplier node ("imposed_T" slot) ──────────────────
    let node0 = nodes[0].id();
    let mut load_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    load_sm.add_cell(&[node0])?;
    load_sm.add_cell(&[mult])?;
    let load_sm = Handle::new(load_sm);
    let mut rhs = SubNodeField::from_poi1(&load_sm, vec!["imposed_T".into(), "q".into()])?;
    rhs.set_value(node0, "q", Q)?;
    rhs.set_value(mult, "imposed_T", T_IMPOSED)?;
    let rhs = NodeField::from_sub(rhs);

    // ── Assemblage + résolution ────────────────────────────────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical solution u(x) = 20 + (Q/k)(1 − x) ─────
    let tol = 1e-10;
    for (i, node) in nodes.iter().enumerate() {
        let x = i as f64 * h;
        let expected = T_IMPOSED + (Q / K) * (1.0 - x);
        let got = solution.value(node.id(), "T")?;
        assert!(
            (got - expected).abs() < tol,
            "T(x={x}) : obtenu {got}, attendu {expected}"
        );
    }
    // La réaction (multiplicateur de Lagrange) équilibre le flux injecté : λ = Q.
    let reaction = solution.value(mult, "lambda_T")?;
    assert!(
        (reaction - Q).abs() < tol,
        "réaction λ : obtenue {reaction}, attendue {Q}"
    );

    Ok(())
}

Exemple Python

La version Python équivalente et documentée est dans le dépôt : examples/thermal_line_1d.py (lancer avec python examples/thermal_line_1d.py après maturin develop). Les compléments 2-D ci-dessous ont eux aussi leur variante Python (thermal_square_2d.py, thermal_convection_2d.py).

Compléments

Exemple : un carré

La généralisation 2-D du cas précédent : un carré \([0,1]^2\) (grille structurée de QUA4), chauffé par une source répartie sur le bord gauche (\(x=0\)) et maintenu à \(T=20\) sur le bord droit (\(x=1\)). Les bords haut et bas ne portent aucune condition : c’est la condition naturelle (flux nul, bord isolé).

Comme les bords latéraux sont isolés, le champ ne dépend pas de \(y\) : le carré redonne le profil de la ligne,

\[ u(x) = 20 + \frac{Q}{k}\,(1 - x), \]

et la réaction totale (somme des multiplicateurs sur le bord imposé) vaut le flux injecté \(Q\).

Mise en donnée d’un flux réparti. Une source répartie se transforme en charges nodales cohérentes \(f_i = \int_\Gamma \varphi\,N_i\,d\Gamma\) par l’opérateur flux — l’analogue de FLUX/PRES de Cast3M. On lui donne le bord (ici un maillage SEG2, intégré comme une ligne : la mesure vient du Jacobien manifold) et la densité de flux (une constante, ou un champ par éléments) ; il renvoie un NodeField sur la composante duale "q", prêt à composer (|) avec le reste du chargement. Sous le capot, pour un flux uniforme sur des éléments linéaires, un nœud intérieur du bord reçoit \(Q\,h\) et un coin \(Q\,h/2\) (somme \(Q\)) — mais on n’a plus à le calculer à la main.

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn thermal_square_recovers_analytical_solution() -> Result<()> {
    // ── Problem data ───────────────────────────────────────────────────────
    const K: f64 = 1.0; // conductivité
    const Q: f64 = 10.0; // flux de chaleur TOTAL injecté sur le bord gauche
    const T_IMPOSED: f64 = 20.0; // température imposée sur le bord droit
    const N: usize = 4; // N×N éléments QUA4
    let h = 1.0 / N as f64;

    // ── Mesh: a structured (N+1)×(N+1) grid of QUA4 on [0,1]² ──────────────
    let coords = Handle::new(Coords::new(2)?);
    let idx = |i: usize, j: usize| j * (N + 1) + i; // nœud colonne i, ligne j
    let mut grid: Vec<Node> = Vec::with_capacity((N + 1) * (N + 1));
    for j in 0..=N {
        for i in 0..=N {
            grid.push(Node::create_in(
                coords.clone(),
                &[i as f64 * h, j as f64 * h],
            )?);
        }
    }
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
    for j in 0..N {
        for i in 0..N {
            mesh.add_cell(&[
                grid[idx(i, j)].id(),
                grid[idx(i + 1, j)].id(),
                grid[idx(i + 1, j + 1)].id(),
                grid[idx(i, j + 1)].id(),
            ])?;
        }
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Dirichlet T = 20 on the right edge (x = 1) ──────────────────────────
    let right_nodes: Vec<Node> = (0..=N).map(|j| grid[idx(N, j)].clone()).collect();
    let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(&right_nodes)?);
    let multiplier = mesh::barycenter(&imposed)?;
    let mults: Vec<Node> = (0..=N)
        .map(|j| multiplier.node(0, j, 0))
        .collect::<Result<_>>()?;

    let conduction = model::heat_conduction(&fes)?;
    let dirichlet = model::dirichlet(&conduction, "T", &imposed, &multiplier, Default::default())?;
    let model = conduction.union(&dirichlet)?;
    // ── Chargement ─────────────────────────────────────────────────────────
    // Source: uniform flux (density Q) on the left edge, turned into
    // consistent nodal loads by the `flux` operator (Cast3m FLUX) — no more
    // Q·h / Q·h/2 distribution by hand. The edge is a SEG2 mesh built on
    // the grid's nodes; it integrates as a line.
    let mut left_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for j in 0..N {
        left_edge.add_cell(&[grid[idx(0, j)].id(), grid[idx(0, j + 1)].id()])?;
    }
    let left_fes = FiniteElementSpace::lagrange1(&left_edge)?;
    let model = model.union(&model::flux(&left_fes, &model, "q".into())?)?;

    let materials =
        pyrucast::ops::element_field::material_field(&model, &[("k", K), ("phi_q", Q)])?;
    let source = pyrucast::ops::node_field::external_forces(&model, &materials)?;

    // Imposed value T = 20 at the multiplier nodes' "imposed_T" slot.
    let mut imposed_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    for m in &mults {
        imposed_sm.add_cell(&[m.id()])?;
    }
    let imposed_sm = Handle::new(imposed_sm);
    let mut imposed_load = SubNodeField::from_poi1(&imposed_sm, vec!["imposed_T".into()])?;
    for m in &mults {
        imposed_load.set_value(m.id(), "imposed_T", T_IMPOSED)?;
    }

    // Loading = the edge's flux + the imposed values (union of the zones).
    let rhs = source.union(&NodeField::from_sub(imposed_load))?;

    // ── Assemblage + résolution ────────────────────────────────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical u(x) = 20 + (Q/k)(1 − x), ∀ y ─────────
    let tol = 1e-9;
    for j in 0..=N {
        for i in 0..=N {
            let x = i as f64 * h;
            let expected = T_IMPOSED + (Q / K) * (1.0 - x);
            let got = solution.value(grid[idx(i, j)].id(), "T")?;
            assert!(
                (got - expected).abs() < tol,
                "T(x={x}, y={}) : obtenu {got}, attendu {expected}",
                j as f64 * h
            );
        }
    }
    // The total reaction on the imposed edge balances the injected flux: Σλ = Q.
    let total_reaction: f64 = mults
        .iter()
        .map(|m| solution.value(m.id(), "lambda_T"))
        .sum::<Result<f64>>()?;
    assert!(
        (total_reaction - Q).abs() < tol,
        "réaction totale : obtenue {total_reaction}, attendue {Q}"
    );

    Ok(())
}

Version Python : examples/thermal_square_2d.py (lancer avec python examples/thermal_square_2d.py après maturin develop).

Conduction orthotrope et anisotrope

Un matériau feuilleté, fibré ou laminé ne conduit pas la chaleur de la même façon dans toutes les directions. La conductivité devient alors un tenseur K, et la rigidité

\[ K_{ij} = \int_\Omega \nabla N_i^{\mathsf T}\, \mathbf{K}\, \nabla N_j \, d\Omega \]

dont le cas isotrope K = k·I redonne le produit scalaire habituel.

C’est le même axe de symétrie matériau qu’en mécanique (chapitre Élasticité orthotrope), avec un tenseur d’ordre 2 au lieu de 4 :

Le repère est donné par des vecteurs — V1X, V1Y en 2-D, V1X…V1Z, V2X…V2Z en 3-D — comme MATE 'DIRECTION' V1 V2 de Cast3M. Ils sont orthonormalisés en interne.

model = pyrucast.model.heat_conduction(fes, symmetry="orthotropic")
materials = pyrucast.element_field.material_field(
    model,
    [("k_1", 12.0), ("k_2", 3.0), ("k_3", 12.0), ("V1X", cos_a), ("V1Y", sin_a)],
)

La conductivité isotrope reste lue au point de Gauss, donc variable à l’intérieur d’une maille ; les constantes orientées sont lues par maille, comme les modules mécaniques.

L’exemple Rust est un test de patch, qui est ce qu’appelle une conductivité orientée : un champ de température linéaire est harmonique pour n’importe quel tenseur constant, donc l’imposer au bord doit le reproduire à l’intérieur quelle que soit K. Le test ne s’arrête pas là — il relit le flux produit et le compare à K·∇T calculé à la main, ce qui est le seul moyen de prendre la rotation en défaut :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::element_field::ElementField;
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::model::Model;
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::symmetry::MaterialSymmetry;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

/// `N×N` QUA4 grid on the unit square.
const N: usize = 3;

#[test]
fn orthotropic_conduction_passes_the_linear_patch_test() -> Result<()> {
    const K1: f64 = 12.0; // along the first material axis
    const K2: f64 = 3.0; // transverse
    let theta = 30.0_f64.to_radians();
    let (c, s) = (theta.cos(), theta.sin());

    let (grid, fes, _coords) = unit_square()?;
    let (model, multipliers) = patch_model(&grid, &fes, MaterialSymmetry::Orthotropic)?;
    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[
            ("k_1", K1),
            ("k_2", K2),
            ("k_3", K1),
            ("V1X", c),
            ("V1Y", s),
        ],
    )?;

    let solution = solve_patch(&model, &materials, &grid, &multipliers)?;

    // The linear field must be reproduced exactly, tensor or no tensor.
    let h = 1.0 / N as f64;
    let tol = 1e-9;
    for j in 0..=N {
        for i in 0..=N {
            let x = i as f64 * h;
            let got = solution.value(grid[j * (N + 1) + i].id(), "T")?;
            assert!((got - x).abs() < tol, "T({x}) = {got}");
        }
    }

    // …and the flux must be the first column of the **rotated** tensor.
    let expect_xx = K1 * c * c + K2 * s * s;
    let expect_yx = (K1 - K2) * c * s;
    let (fx, fy) = uniform_flux(&model, &solution, &fes, &materials, &grid)?;
    assert!(
        (fx - expect_xx).abs() < 1e-9,
        "flux_x = {fx}, expected {expect_xx}"
    );
    assert!(
        (fy - expect_yx).abs() < 1e-9,
        "flux_y = {fy}, expected {expect_yx}"
    );
    Ok(())
}

Avec ∇T = (1, 0), le flux est la première colonne de K : K_xx = k₁cos²θ + k₂sin²θ et K_yx = (k₁ − k₂)·cosθ·sinθ. Le terme extra-diagonal n’est non nul que si le matériau est à la fois anisotrope et désaligné — précisément le cas qu’une rotation fausse manquerait.

Rayonnement à l’infini (Stefan-Boltzmann)

Une surface qui échange avec un environnement lointain à \(T_\infty\) rayonne

\[ q\cdot n = \sigma\,\varepsilon\,\big(T^4 - T_\infty^4\big) \]

où \(\sigma\) est la constante de Stefan-Boltzmann et \(\varepsilon\) l’émissivité. Primale "T", duale "q" — les mêmes degrés de liberté que la conduction, donc un bord rayonnant se couple directement dans sa rigidité, comme la convection. Et comme elle, il n’a besoin d’aucune normale : la direction est déjà consommée en écrivant q·n, il ne reste sous l’intégrale qu’un scalaire et la mesure de surface.

Ce qui change par rapport à la convection : c’est non linéaire

La loi de Newton est linéaire en T, si bien que la convection ne contribue qu’une matrice de film constante. T⁴ ne l’est pas, d’où trois termes :

Linéariser la rigidité autour de \(T_\infty\) plutôt qu’autour de l’état courant est ce qui la laisse être une matrice constante : c’est l’opérateur dont on part pour une boucle de Newton, et à lui seul une itération de Picard tout à fait utilisable. La tangente porte la vraie non-linéarité : elle évalue 4σεT³ à la température courante, au point de Gauss, quand on la demande — comme le D_alg plastique, et pour la même raison : personne d’autre ne la lirait.

Deux natures

Le rayonnement déclare [Thermal, Radiation]. Un bord rayonnant fait partie du problème thermique — filter("thermal") doit le rendre — tandis que filter("radiation") isole le terme non linéaire à part, pour l’assembler ou l’inspecter seul. C’est le premier usage du caractère ensembliste de physics().

Unités

sigma vaut par défaut la constante SI, et T est alors une température absolue (Kelvin) : une puissance quatrième n’a aucune invariance permettant de translater une origine. Dans un autre système d’unités, fournir sigma comme composante matériau.

conduction = pyrucast.model.heat_conduction(volume)
model = conduction | pyrucast.model.radiation(bord, conduction)
materials = pyrucast.element_field.material_field(
    model, [("k", 20.0), ("emis", 0.8), ("T_inf", 300.0)]
)

Ce que ça vaut comme vérification

Deux choses se contrôlent sans acrobatie analytique : le flux rayonné doit valoir exactement σε(T⁴ − T_∞⁴) fois l’aire, et la tangente doit être la dérivée du résidu. Une loi en T⁴ est précisément là où une tangente incohérente se cache — Newton ramperait au lieu de converger quadratiquement — d’où sa comparaison à une différence finie :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::element_field::ElementField;
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::model::Model;
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::radiation::STEFAN_BOLTZMANN;
use pyrucast::models::Physics;
use pyrucast::ops::element_field;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

const EMIS: f64 = 0.8; // emissivity
const T_INF: f64 = 300.0; // far-field temperature (K)
const T_WALL: f64 = 500.0; // temperature imposed on the far side (K)

#[test]
fn the_radiated_flux_matches_stefan_boltzmann() -> Result<()> {
    let (fixture, materials) = radiating_square()?;

    // The boundary is at the uniform wall temperature: interpolate it to the
    // Gauss points, integrate the law, and scatter back to the nodes.
    let temperature = uniform_temperature(&fixture, T_WALL)?;
    let at_gauss = element_field::interp_to_gauss(&temperature, &fixture.boundary_fes)?;
    let state =
        element_field::behavior::integrate(&fixture.radiation, &at_gauss, None, &materials, None)?;
    let reaction = pyrucast::ops::node_field::internal_forces(
        &fixture.radiation,
        &state,
        &temperature,
        &materials,
    )?;

    // The radiating edge has unit length, so the total flux is the density.
    let expected = STEFAN_BOLTZMANN * EMIS * (T_WALL.powi(4) - T_INF.powi(4));
    let total: f64 = fixture
        .edge
        .iter()
        .map(|n| reaction.value(n.id(), "q").unwrap_or(0.0))
        .sum();
    assert!(
        (total - expected).abs() < 1e-9 * expected.abs(),
        "radiated flux {total}, expected {expected}"
    );
    Ok(())
}

Convection de surface (Robin / film)

Le modèle BoundaryTransfer (src/models/boundary_transfer.rs) ajoute un échange convectif avec un fluide à température ambiante \(T_\text{ext}\) sur un bord : la loi de Newton du refroidissement

\[ q\cdot n = h\,\big(T - T_\text{ext}\big) \]

où \(h\) est le coefficient d’échange (film). Injectée dans le terme de bord de la forme faible de la conduction, elle se scinde en deux ingrédients :

\[ \underbrace{K_{ij} = h \int_\Gamma N_i\,N_j\,d\Gamma}{\text{matrice de film (raideur)}} \qquad \underbrace{f_i = h\,T\text{ext} \int_\Gamma N_i\,d\Gamma}_{\text{charge (second membre)}} \]

On le construit contre la conduction qu’il refroidit, en lui passant les couples de variables à échanger — ceux de la conduction, ce qui fait que le terme se couple directement dans sa raideur. La conduction, elle, lui donne sa nature thermique, et refuse un couple qu’elle n’assemble pas :

conduction = pyrucast.model.heat_conduction(bord_fes)
film = pyrucast.model.boundary_transfer(bord_fes, conduction, [("T", "q")])

Ce modèle n’a rien de thermique : la même loi décrit un transfert de masse en surface ou une fondation élastique, selon les composantes qu’on lui donne, et il partage son noyau avec le transfert d’interface. Voir Échanges pour la loi commune, la structure en quatre blocs et le choix entre un échange et une contrainte.

Mise en donnée. Le modèle fournit la matrice de film ; la part externe \(h\,T_\text{ext}\) est un chargement, bâti avec le même opérateur flux que la source (densité \(h\,T_\text{ext}\)). Le terme de film rend la matrice définie : un problème purement Neumann + convection est bien posé sans Dirichlet.

Exemple. Une dalle \([0,1]^2\) chauffée par un flux \(Q\) sur le bord gauche et refroidie par convection sur le bord droit (haut/bas isolés). Tout le flux ressort par convection, d’où le profil linéaire

\[ T(x) = T_\text{ext} + \frac{Q}{h} + \frac{Q}{k}\,(1 - x). \]

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn thermal_convection_recovers_analytical_solution() -> Result<()> {
    // ── Problem data ───────────────────────────────────────────────────────
    const K: f64 = 2.0; // conductivité
    const Q: f64 = 10.0; // densité de flux injectée sur le bord gauche
    const H: f64 = 5.0; // coefficient d'échange (film) sur le bord droit
    const T_EXT: f64 = 20.0; // température ambiante du fluide
    const N: usize = 4; // N×N éléments QUA4
    let step = 1.0 / N as f64;

    // ── Mesh: a structured (N+1)×(N+1) grid of QUA4 on [0,1]² ──────────────
    let coords = Handle::new(Coords::new(2)?);
    let idx = |i: usize, j: usize| j * (N + 1) + i; // nœud colonne i, ligne j
    let mut grid: Vec<Node> = Vec::with_capacity((N + 1) * (N + 1));
    for j in 0..=N {
        for i in 0..=N {
            grid.push(Node::create_in(
                coords.clone(),
                &[i as f64 * step, j as f64 * step],
            )?);
        }
    }
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
    for j in 0..N {
        for i in 0..N {
            mesh.add_cell(&[
                grid[idx(i, j)].id(),
                grid[idx(i + 1, j)].id(),
                grid[idx(i + 1, j + 1)].id(),
                grid[idx(i, j + 1)].id(),
            ])?;
        }
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : conduction (volume) + convection (bord droit x = 1) ───────
    // Le bord droit est un maillage SEG2 bâti sur les nœuds de la grille ;
    // it integrates as a line (film matrix h ∫ N_i N_j dΓ).
    let mut right_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for j in 0..N {
        right_edge.add_cell(&[grid[idx(N, j)].id(), grid[idx(N, j + 1)].id()])?;
    }
    let right_fes = FiniteElementSpace::lagrange1(&right_edge)?;

    let conduction = model::heat_conduction(&fes)?;
    let convection =
        model::boundary_transfer(&right_fes, &conduction, vec![("T".into(), "q".into())])?;
    let model = conduction.union(&convection)?;

    // Material: k for the conduction, h and the ambient for the convection (each
    // sub-model takes the component it requires from the supplied list).
    // ── Chargement ─────────────────────────────────────────────────────────
    // Source: uniform flux (density Q) on the left edge, as nodal loads
    // cohérentes via `flux`.
    let mut left_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for j in 0..N {
        left_edge.add_cell(&[grid[idx(0, j)].id(), grid[idx(0, j + 1)].id()])?;
    }
    let left_fes = FiniteElementSpace::lagrange1(&left_edge)?;
    let model = model.union(&model::flux(&left_fes, &model, "q".into())?)?;

    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[("k", K), ("h_T", H), ("a_ext_T", T_EXT), ("phi_q", Q)],
    )?;

    // Both given terms — the left edge's source and the convection's external
    // part h·T_ext — belong to the model, which returns them together. Nothing
    // left to union by hand, hence nothing left to forget.
    let rhs = pyrucast::ops::node_field::external_forces(&model, &materials)?;

    // ── Assembly + solve (K made definite by the film term) ────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical T(x) = T_ext + Q/h + (Q/k)(1 − x), ∀ y ─
    let tol = 1e-9;
    for j in 0..=N {
        for i in 0..=N {
            let x = i as f64 * step;
            let expected = T_EXT + Q / H + (Q / K) * (1.0 - x);
            let got = solution.value(grid[idx(i, j)].id(), "T")?;
            assert!(
                (got - expected).abs() < tol,
                "T(x={x}, y={}) : obtenu {got}, attendu {expected}",
                j as f64 * step
            );
        }
    }

    // Energy balance: all the injected flux leaves by convection, so the right
    // edge's temperature is exactly T_ext + Q/h.
    let t_right = solution.value(grid[idx(N, 0)].id(), "T")?;
    assert!(
        (t_right - (T_EXT + Q / H)).abs() < tol,
        "T(x=1) : obtenu {t_right}, attendu {}",
        T_EXT + Q / H
    );

    Ok(())
}

Version Python : examples/thermal_convection_2d.py (lancer avec python examples/thermal_convection_2d.py après maturin develop).