Rayonnement

En bref

Le rayonnement d’une surface vers un environnement lointain à la température \(T_\infty\) (loi de Stefan-Boltzmann). Il se pose sur le bord d’une conduction et partage ses variables : il entre directement dans sa matrice.

fespace est l’espace du bord, target la conduction qu’il refroidit. Une cible qui n’assemble pas T avec sa duale q est refusée : le rayonnement écrirait dans une ligne que personne ne résout.

Les deux tags permettent de garder le rayonnement dans le problème thermique (tag thermal) ou de le prendre seul (tag radiation) pour l’assembler ou l’inspecter à part.

Théorie

Le flux qui quitte la surface vaut

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

où \(\sigma\) est la constante de Stefan-Boltzmann et \(\varepsilon\) l’émissivité. Aucune normale n’entre dans le calcul : la direction est déjà dans \(q\cdot n\), il ne reste qu’un scalaire et la mesure de la surface.

La loi n’est pas linéaire en \(T\), d’où trois termes :

Linéarisée autour de \(T_\infty\), la raideur est une matrice constante : c’est l’opérateur de départ d’une boucle de Newton, et à elle seule un opérateur de point fixe utilisable. La tangente porte la non-linéarité : elle est évaluée à la température de chaque point de Gauss.

sigma vaut par défaut la constante SI, 5,670374419·10⁻⁸ W·m⁻²·K⁻⁴, et T est alors une température absolue, en kelvins : une puissance quatrième ne permet pas de déplacer l’origine. Dans un autre système d’unités, on donne sigma dans le matériau.

Variables

Ce qu’il produit

Dans un calcul

Le flux se calcule aux points de Gauss de la face. La température y est interpolée à partir de ses valeurs aux nœuds, \( T_g = \sum_j N_j(x_g)\,T_j \), élevée à la puissance quatrième en ce point, et le flux \( \sigma\varepsilon(T_g^4 - T_\infty^4) \) est intégré avec le poids \( w_g \) de chaque point (la mesure de la surface comprise) :

\[ R_i = \sum_g w_g\, N_i(x_g)\, \sigma\varepsilon\,(T_g^4 - T_\infty^4), \qquad (K_r)_{ij} = \sum_g w_g\, 4\sigma\varepsilon\,T_g^3\, N_i(x_g)\, N_j(x_g) \]

\( R \) est le résidu, que rendent les forces internes, et \( K_r \) sa dérivée exacte par rapport aux températures nodales : la tangente.

Le calcul pas à pas annule le résidu par corrections successives. À chaque pas de temps, pyrucast.stepping part de la température du pas précédent et calcule le résidu de tout le problème — conduction, capacité, chargements, rayonnement. Il résout la correction avec un opérateur qui contient la tangente \( K_r \), prise une fois au début du pas, puis recommence avec le résidu de la nouvelle température, jusqu’à ce que la correction tombe sous precision (\( 10^{-4} \) par défaut) relativement à la température. Le résidu suit donc la température à chaque itération ; l’opérateur, lui, reste celui du début du pas. Comme il contient la dérivée de \( T^4 \), une ou deux corrections suffisent quand le pas de temps est raisonnable (algorithme détaillé dans Calcul pas à pas, Algo 2).

La raideur, film linéarisé autour de \( T_\infty \), est une matrice constante : employée seule à chaque itération, elle donne une itération de point fixe, qui converge plus lentement que la tangente.

Une autre méthode consiste à écrire le rayonnement comme une convection dont le coefficient dépend de la température. Puisque \( T^4 - T_\infty^4 = (T^2 + T_\infty^2)(T + T_\infty)(T - T_\infty) \), le flux vaut \( h_r\,(T - T_\infty) \), avec

\[ h_r = \sigma\varepsilon\,(T^2 + T_\infty^2)(T + T_\infty). \]

On calcule \( h_r \) aux nœuds, avec la température nodale de l’itération en cours, on l’interpole aux points de Gauss, et chaque itération résout un problème de convection : \( \int_\Gamma h_r N_i N_j \) dans la matrice, \( \int_\Gamma h_r\,T_\infty N_i \) au second membre. On recommence avec la nouvelle température jusqu’à ce qu’elle ne bouge plus.

Deux choses changent. \( h_r \) est une sécante, pas la tangente \( 4\sigma\varepsilon T^3 \) : chaque itération corrige moins, et il en faut davantage pour la même précision. Et le problème discret n’est pas tout à fait le même : \( h_r \) interpolé depuis les nœuds n’est pas \( h_r \) calculé avec la température des points de Gauss. Les deux méthodes convergent vers la même solution quand le maillage s’affine ; sur un maillage donné, leurs solutions convergées diffèrent de l’ordre de l’erreur de discrétisation.

Exemple

Une plaque à 1000 K qui rayonne par son bord supérieur vers 300 K. On lit le flux au bord, puis l’écart entre tangente et raideur au coin du bord : un seul segment, de longueur L = 0,5, touche ce coin, où \(\int N_0^2 = L/3\) :

conduction = pyrucast.model.heat_conduction(plaque)
rayonnement = pyrucast.model.radiation(pyrucast.FiniteElementSpace(haut), conduction)
print(rayonnement[0].physics())  # ['thermal', 'radiation']
modele = conduction | rayonnement
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 20.0), ("emis", 0.8), ("T_inf", 300.0)]
)

T = pyrucast.NodeField(carre, ["T"])
T.unit().set_uniform("T", 1000.0)  # en kelvins
entree = pyrucast.element_field.strain(modele, T, materiaux)  # ∇T dedans, T au bord
repos = pyrucast.element_field.initial_state(modele, materiaux)
sortie = pyrucast.element_field.integrate_behavior(modele, entree, materiaux, repos)
print(sortie[1].value(0, 0, "flux"))  # 44995.55… = σε(T⁴ − T∞⁴)

K = pyrucast.matrix.stiffness(modele, materiaux)  # film linéarisé en T∞
Kt = pyrucast.matrix.tangent(
    modele, materiaux, entree, repos
)  # film à la température T
coin = haut.node(0, 0, 0)
print(
    Kt.get(coin, "q", coin, "T") - K.get(coin, "q", coin, "T")
)  # 29.425… = 4σε(T³ − T∞³)·L/3

Côté Rust

Radiation (src/models/radiation.rs) garde son sous-espace de bord. Son noyau de raideur évalue 4σεT∞³ au point de Gauss, son comportement (integrate_point) le flux, et tangent_point la dérivée 4σεT³. STEFAN_BOLTZMANN est la constante par défaut.

Le test compare le flux rayonné à σε(T⁴ − T∞⁴) fois l’aire, et la tangente à une différence finie du résidu :

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::symmetry::MaterialSymmetry;
use pyrucast::models::Physics;
use pyrucast::ops::element_field;
use pyrucast::ops::model;
use pyrucast::ops::solver::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,
        &element_field::behavior::initial_state(&fixture.radiation, &materials)?,
        &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(())
}