Échanges et flux imposé

Trois physiques sans comportement, qui s’appliquent à une autre physique, leur cible, et en prennent les variables et le tag :

• l’échange de bord, entre un bord et un milieu extérieur connu ;
• le transfert d’interface, entre deux maillages qui se font face ;
• le flux imposé, une densité donnée sur un bord ou un volume.

Les deux échanges suivent la même loi : un flux proportionnel à un écart,

\[ q\cdot n = h\,(a - b). \]

Pour l’échange de bord, \(b\) est une donnée, la valeur ambiante ; pour le transfert d’interface, c’est une inconnue, le champ de l’autre maillage. Les deux partagent leur noyau, transfer::exchange_matrix.

Ce qui est échangé est donné par l’appelant. On passe la liste des variables échangées (["T"], ["c_H2"], ["u_x", "u_y"]). Chacune doit être une inconnue de la cible, qui lui donne sa duale, la ligne où le flux est écrit. Le coefficient \(h\) est un par variable, nommé d’après elle : h_T, h_c_H2, h_u_x. L’échange prend le tag de la cible : thermal sur une conduction, diffusion sur une diffusion, mechanical sur une élasticité. Un nom inconnu de la cible, ou des variables de deux tags dans un même échange, sont refusés à la construction.

Échange de bord

En bref

Un échange entre un bord et un milieu extérieur : la convection avec un fluide (loi de Newton), le transfert de masse en surface, ou, sur les déplacements, une fondation élastique.

Théorie

Sur le bord \(\Gamma\), le flux sortant vaut \(h\,(a - a_\text{ext})\). Dans la forme faible, il donne une matrice de film et un second membre :

\[ K_{ij} = \int_\Gamma h\,N_i\,N_j\,d\Gamma , \qquad f_i = \int_\Gamma h\,a_\text{ext}\,N_i\,d\Gamma . \]

Aucune normale n’entre dans le calcul : elle est déjà dans \(q\cdot n\), et il ne reste sous l’intégrale qu’un scalaire et la mesure du bord, qui ne dépend pas de l’orientation des mailles. Le film rend la matrice définie : un problème avec un flux imposé et un échange, sans aucune température imposée, a une solution unique.

Variables

Pour la variable T d’une conduction :

Ce qu’il produit

Exemple

Un carré de conductivité 2, avec une source de 10 sur le bord gauche (un flux imposé) et un film h = 5 vers 20 sur le bord droit. Tout ce qui entre ressort par le film, d’où \(T = T_\text{ext} + Q/h\) à droite et \(T_\text{ext} + Q/h + Q/k\) à gauche :

import pyrucast

# Un carré [0, 1]² : ses bords gauche et droit sont les deux lignes balayées.
c = pyrucast.Coords(2)
gauche = pyrucast.mesh.line(c.add_node([0.0, 0.0]), c.add_node([0.0, 1.0]), 2)
droite = pyrucast.mesh.line(c.add_node([1.0, 0.0]), c.add_node([1.0, 1.0]), 2)
carre = pyrucast.mesh.sweep(gauche, droite, 2)

conduction = pyrucast.model.heat_conduction(pyrucast.FiniteElementSpace(carre))
source = pyrucast.model.flux(pyrucast.FiniteElementSpace(gauche), conduction, "q")
film = pyrucast.model.boundary_transfer(
    pyrucast.FiniteElementSpace(droite), conduction, ["T"]
)
print(source[0].physics(), film[0].physics())  # ['thermal'] ['thermal']
modele = conduction | source | film
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 2.0), ("phi_q", 10.0), ("h_T", 5.0), ("a_ext_T", 20.0)]
)

second_membre = pyrucast.node_field.external_forces(
    modele, materiaux
)  # φ à gauche, h·T_ext à droite
T = pyrucast.matrix.stiffness(modele, materiaux).solve(second_membre)
print(round(T.value(droite.node(0, 0, 0), "T"), 9))  # 22.0 = T_ext + Q/h
print(round(T.value(gauche.node(0, 0, 0), "T"), 9))  # 27.0 = T_ext + Q/h + Q/k

Le même échange sur les déplacements est une fondation élastique : un ressort réparti par direction.

plaque = pyrucast.model.elasticity(pyrucast.FiniteElementSpace(carre), "plane_stress")
semelle = pyrucast.model.boundary_transfer(
    pyrucast.FiniteElementSpace(gauche), plaque, ["u_x", "u_y"]
)
print(semelle[0].physics())  # ['mechanical']
print(semelle[0].material_components())  # ['h_u_x', 'h_u_y', 'a_ext_u_x', 'a_ext_u_y']

Transfert d’interface

En bref

Un échange entre deux maillages qui se font face sans partager leurs nœuds : une résistance de contact, un revêtement, un joint, une membrane. Le champ saute à la traversée, et le flux qui passe est proportionnel au saut.

side_a et side_b sont les espaces des deux faces. Elles doivent être conformes : même type de maille, même nombre de mailles, la maille i face à la maille i, et leur nœud local k à moins de tol l’un de l’autre. C’est vérifié à la construction, et une interface non conforme est refusée. Un sous-modèle est créé par paire de zones des deux côtés, prises dans l’ordre.

Théorie

La forme faible du terme d’interface est

\[ \int_\Gamma h\,(a_1 - a_2)\,(\delta a_1 - \delta a_2)\,d\Gamma , \]

qui donne quatre blocs sur les inconnues des deux côtés :

\[ \begin{bmatrix} +K & -K \\ -K & +K \end{bmatrix}, \qquad K_{ij} = \int_\Gamma h\,N_i\,N_j\,d\Gamma . \]

Les deux blocs diagonaux sont deux films, un par côté. Les deux autres ont leurs lignes sur un maillage et leurs colonnes sur l’autre ; chacun seul n’est pas symétrique, leur paire l’est. Les variables échangées ne se couplent pas entre elles : la chaleur qui traverse un joint ne pousse pas l’hydrogène.

C’est la version par pénalité de ce qu’un multi-points impose exactement : quand \(h\) grandit, le saut tend vers zéro. Si \(h\) est une donnée physique (une résistance de contact mesurée), c’est un transfert ; s’il n’est là que pour lier deux surfaces, c’est une contrainte.

Variables

Pour la variable T :

Ce qu’il produit

Exemple

Deux carrés côte à côte, un joint h = 4 entre eux, un flux de 1 entrant à gauche et T = 0 à droite. Le saut à la traversée vaut \(Q/h\) :

# Deux carrés côte à côte ; la ligne x = 1 existe deux fois, avec des nœuds distincts.
c = pyrucast.Coords(2)
x0 = pyrucast.mesh.line(c.add_node([0.0, 0.0]), c.add_node([0.0, 1.0]), 2)
x1_g = pyrucast.mesh.line(c.add_node([1.0, 0.0]), c.add_node([1.0, 1.0]), 2)
x1_d = pyrucast.mesh.line(c.add_node([1.0, 0.0]), c.add_node([1.0, 1.0]), 2)
x2 = pyrucast.mesh.line(c.add_node([2.0, 0.0]), c.add_node([2.0, 1.0]), 2)
corps = pyrucast.model.heat_conduction(
    pyrucast.FiniteElementSpace(pyrucast.mesh.sweep(x0, x1_g, 2))
) | pyrucast.model.heat_conduction(
    pyrucast.FiniteElementSpace(pyrucast.mesh.sweep(x1_d, x2, 2))
)
joint = pyrucast.model.interface_transfer(
    pyrucast.FiniteElementSpace(x1_g), pyrucast.FiniteElementSpace(x1_d), corps, ["T"]
)
source = pyrucast.model.flux(pyrucast.FiniteElementSpace(x0), corps, "q")
fixe = x2.to_poi1()  # les trois nœuds de x = 2
appui = pyrucast.model.dirichlet(fixe, corps, "T")
modele = corps | joint | source | appui
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 1.0), ("h_T", 4.0), ("phi_q", 1.0)]
)

second_membre = pyrucast.node_field.external_forces(
    modele, materiaux
) | appui.constraint_rhs([(fixe.node(0, i, 0), 0.0) for i in range(3)])
T = pyrucast.matrix.stiffness(modele, materiaux).solve(second_membre)
saut = T.value(x1_g.node(0, 0, 0), "T") - T.value(x1_d.node(0, 0, 0), "T")
print(round(saut, 9))  # 0.25 = Q/h : le champ saute à la traversée

Flux imposé

En bref

Une densité donnée, intégrée en forces nodales : une source de chaleur sur un bord ou dans un volume, une traction, une pression. Elle ne dépend pas de la solution, et n’entre dans aucune matrice.

dual est la ligne de la cible où la densité est écrite : "q" pour une source de chaleur, "f_x" pour une traction en x. "n" désigne une pression : la densité agit selon la normale de la maille et s’écrit dans toutes les lignes de déplacement (f_x, f_y, f_z). La normale est celle de la configuration active : sur une configuration déformée, la pression suit la paroi.

Théorie

\[ f_i = \int \varphi\,N_i\,d\Gamma , \qquad f_{i,a} = \int \varphi\,n_a\,N_i\,d\Gamma \;\text{pour une pression.} \]

Variables

La densité est une composante du champ matériau : uniforme, ou variable d’un point de Gauss à l’autre.

Ce qu’il produit

Exemple

La source de l’échange de bord ci-dessus est un flux imposé : model.flux(gauche, conduction, "q") et la composante phi_q du matériau.

Côté Rust

BoundaryTransfer (src/models/boundary_transfer.rs) et InterfaceTransfer (src/models/interface_transfer.rs) s’appuient sur src/models/transfer.rs, qui porte leur noyau commun, exchange_matrix, la vérification des variables contre la cible et le choix du tag (owner_physics). Ils implémentent Domain sans Behavior. Ce qu’ils lisent pour leurs forces internes est la primale interpolée aux points de Gauss (residual_input), et non un état.

Le transfert d’interface rend ses quatre blocs par contributions : deux blocs ordinaires et deux blocs de couplage (Contribution::Coupling), dont les lignes et les colonnes vivent sur deux sous-espaces différents.

Flux (src/models/flux.rs) ne contribue à aucune matrice ; il déclare seulement son terme de forces externes (external_force_contribution).

Le test de la fondation élastique compare le déplacement d’une barre poussée sur un ressort réparti à \(u = q / (E/L + h)\), sur quatre décades de raideur :

/// A boundary exchange on the displacement is an **elastic foundation**: the
/// face rests on a distributed spring, and the analytical `u = q/(E/L + h)`
/// comes out for stiffnesses spanning four decades.
#[test]
fn a_boundary_exchange_on_displacements_is_an_elastic_foundation() -> Result<()> {
    for h in [0.0, 1e3, 1e4, 1e5, 1e6] {
        let u = free_face_displacement(4, h)?;
        let exact = Q / (E / L + h);
        assert!(
            (u - exact).abs() < 1e-9 * exact.abs().max(1e-12),
            "h = {h}: face displacement {u}, exact {exact}"
        );
    }
    Ok(())
}

Le test du transfert d’interface vérifie le profil de part et d’autre et le saut \(q/h\) :

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::interface_transfer::DEFAULT_TOL;
use pyrucast::models::symmetry::MaterialSymmetry;
use pyrucast::models::{Physics, IMPOSED};
use pyrucast::ops::model;
use pyrucast::ops::solver::solve;
use pyrucast::Result;

/// The diffusing species — every name of the Fick physics carries it.
const SPECIES: &str = "H2";

const D: f64 = 2.0; // diffusivity, both squares
const Q: f64 = 10.0; // flux density injected at x = 0
const C_RIGHT: f64 = 1.0; // concentration imposed at x = 2

#[test]
fn an_interface_law_makes_the_field_jump() -> Result<()> {
    const H: f64 = 5.0; // transfer coefficient
    let (geom, solution) = solve_two_squares(H)?;

    let c = |n: &Node| solution.value(n.id(), &format!("c_{SPECIES}"));
    let far_left = c(&geom.left[0])?; // (0, 0)
    let left_face = c(&geom.left[1])?; // (1, 0), left side
    let right_face = c(&geom.right[0])?; // (1, 0), right side
    let far_right = c(&geom.right[1])?; // (2, 0)

    let tol = 1e-10;
    assert!((far_right - C_RIGHT).abs() < tol, "c(2) = {far_right}");
    // Slope q/D over the unit width of each square.
    assert!(
        (far_left - left_face - Q / D).abs() < tol,
        "left square: {far_left} → {left_face}"
    );
    assert!(
        (right_face - far_right - Q / D).abs() < tol,
        "right square: {right_face} → {far_right}"
    );
    // …and the jump across the interface is q/h — the exchange law itself.
    let jump = left_face - right_face;
    assert!(
        (jump - Q / H).abs() < tol,
        "jump = {jump}, expected {}",
        Q / H
    );
    Ok(())
}