Conduction thermique

En bref

La conduction de la chaleur dans un solide : la température T est l’inconnue, et la chaleur circule des zones chaudes vers les zones froides.

symmetry choisit la forme de la conductivité : "isotropic" (un scalaire k), "orthotropic" ou "anisotropic" (un tenseur orienté par des axes matériau).

Théorie

En régime transitoire, l’équation de la chaleur s’écrit

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

avec \(\mathbf K\) le tenseur de conductivité, \(\rho\,c_p\) la capacité volumique et \(Q\) une source. En régime permanent, le premier terme disparaît. Multipliée par une température virtuelle \(v\) et intégrée par parties, elle devient

\[ \int_\Omega \nabla v\cdot\mathbf K\,\nabla T\,d\Omega + \int_\Omega v\,\rho\,c_p\,\dot T\,d\Omega = \int_\Omega v\,Q\,d\Omega + \int_\Gamma v\,\bar q\,d\Gamma , \]

où \(\bar q\) est le flux qui entre par le bord. Un bord sans condition est isolé : rien n’y entre.

Avec \(T \approx \sum_j N_j\,T_j\), chaque maille donne une matrice de conduction et une matrice de capacité :

\[ K_{ij} = \int_e \nabla N_i^{\mathsf T}\,\mathbf K\,\nabla N_j\,d\Omega , \qquad C_{ij} = \int_e \rho\,c_p\,N_i\,N_j\,d\Omega , \]

intégrées aux points de Gauss, où le matériau est lu : la conductivité et la capacité peuvent varier à l’intérieur d’une maille. Le système semi-discret est \(C\,\dot T + K\,T = F\) ; un pas d’Euler implicite résout \((C/\Delta t + K)\,T_{n+1} = F + (C/\Delta t)\,T_n\).

Conductivité orientée. Un matériau fibré ou feuilleté conduit mieux dans certaines directions. Le tenseur \(\mathbf K\) est alors construit dans les axes matériau, diagonal (orthotropic, conductivités principales k_1, k_2, k_3) ou complet (anisotropic, k_11…k_33), puis tourné dans les axes globaux. Les axes sont des vecteurs du champ matériau : V1X, V1Y en 2-D, le second axe étant la normale dans le plan ; V1X…V1Z et V2X…V2Z en 3-D, le troisième étant V1 × V2. Ils sont normés en interne et peuvent varier d’une maille à l’autre. En 2-D, k_3 (la direction normale au plan) est demandé aussi : le tenseur est construit en 3-D, puis réduit au plan.

Variables

La sortie est 𝐊·∇T, le flux de la forme faible. Le flux de Fourier, celui qui va du chaud vers le froid, est son opposé : −𝐊·∇T.

Ce qu’elle produit

Exemple

Une barre de conductivité 2, chauffée par un flux de 10 à gauche, à 20 à droite. Le profil est linéaire, \(T(x) = 20 + Q\,(1 - x)/k\) :

import pyrucast

# Une barre [0, 1] de quatre SEG2 : un flux Q = 10 entre à gauche, T = 20 à droite.
c = pyrucast.Coords(1)
gauche = c.add_node([0.0])
droite = c.add_node([1.0])
barre = pyrucast.mesh.line(gauche, droite, 4)
fes = pyrucast.FiniteElementSpace(barre)

conduction = pyrucast.model.heat_conduction(fes)
appui = pyrucast.model.dirichlet(
    pyrucast.mesh.poi1_from_nodes([droite]), conduction, "T"
)
modele = conduction | appui
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 2.0), ("rho", 1.0), ("cp", 3.0)]
)

source = pyrucast.NodeField(barre, ["q"])  # sur les nœuds de la barre
source.unit().set_value(gauche, "q", 10.0)
second_membre = source | appui.constraint_rhs([(droite, 20.0)])

K = pyrucast.matrix.stiffness(modele, materiaux)
T = K.solve(second_membre)
print(T.value(gauche, "T"))  # 25.0 = 20 + Q·L/k
print(T.value(droite, "T"))  # 20.0

Le gradient et le flux aux points de Gauss, et les forces internes, qui rendent la source :

gradient = pyrucast.element_field.strain(modele, T, materiaux)
print(gradient.components())  # ['grad_T_x']
print(gradient[0].value(0, 0, "grad_T_x"))  # -5.0

repos = pyrucast.element_field.initial_state(
    modele, materiaux
)  # l'état de départ de la loi
flux = pyrucast.element_field.integrate_behavior(modele, gradient, materiaux, repos)
print(flux[0].value(0, 0, "flux_x"))  # -10.0 = k·∂T/∂x

f_int = pyrucast.node_field.internal_forces(modele, flux, T, materiaux)
print(f_int.value(gauche, "q"))  # 10.0 : la source que la conduction équilibre

Un pas de temps, depuis une barre à 20 partout :

C = pyrucast.matrix.mass(modele, materiaux)  # la capacité ρ·cp ; l'appui n'y entre pas
T0 = pyrucast.NodeField(barre, ["T"])
T0.unit().set_uniform("T", 20.0)

dt = 0.01
A = C / dt + K  # Euler implicite : (C/Δt + K)·T₁ = F + (C/Δt)·T₀
T1 = A.solve(second_membre + (C / dt) * T0)
print(round(T1.value(gauche, "T"), 3))  # 20.306 : la chaleur commence à entrer

Une conductivité orthotrope, d’axe principal (0,6 ; 0,8) :

c2 = pyrucast.Coords(2)
bas = pyrucast.mesh.line(c2.add_node([0.0, 0.0]), c2.add_node([1.0, 0.0]), 2)
haut = pyrucast.mesh.line(c2.add_node([0.0, 1.0]), c2.add_node([1.0, 1.0]), 2)
carre = pyrucast.mesh.sweep(bas, haut, 2)
plaque = pyrucast.FiniteElementSpace(carre)

orientee = pyrucast.model.heat_conduction(plaque, symmetry="orthotropic")
print(orientee[0].material_components())  # ['k_1', 'k_2', 'k_3', 'V1X', 'V1Y']
materiaux = pyrucast.element_field.material_field(
    orientee,
    [("k_1", 12.0), ("k_2", 3.0), ("k_3", 3.0), ("V1X", 0.6), ("V1Y", 0.8)],
)
K = pyrucast.matrix.stiffness(orientee, materiaux)

Côté Rust

HeatConduction (src/models/heat_conduction.rs) garde son sous-espace et sa MaterialSymmetry. Les matrices élémentaires sont des fonctions pures, element_stiffness et element_capacity, qui lisent le matériau ligne par ligne aux points de Gauss par les indices résolus une fois par zone (ElementLayout). Le tenseur orienté est construit par src/models/symmetry.rs, le même code que pour la diffusion et l’élasticité.

Le même problème que l’exemple, en Rust :

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::models::symmetry::MaterialSymmetry;
use pyrucast::ops::model;
use pyrucast::ops::solver::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 constraint mints its multiplier beside the imposed node (a fresh
    // co-located node); its value is written there.
    let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(
        nodes.last().unwrap(),
    ))?);

    let conduction = model::heat_conduction(&fes, MaterialSymmetry::Isotropic)?;
    let dirichlet = model::dirichlet(&imposed, &conduction, "T", Default::default())?;
    let mult = dirichlet.multiplier_mesh()?.node(0, 0, 0)?.id();
    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" 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".into(), "q".into()])?;
    rhs.set_value(node0, "q", Q)?;
    rhs.set_value(mult, "imposed", 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")?;
    assert!(
        (reaction - Q).abs() < tol,
        "réaction λ : obtenue {reaction}, attendue {Q}"
    );

    Ok(())
}

Le test de la conduction orientée impose un champ linéaire au bord, qui doit être reproduit à l’intérieur pour tout tenseur constant, puis relit le flux et le compare à 𝐊·∇T calculé à la main :

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::model;
use pyrucast::ops::solver::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(())
}