Diffusion (loi de Fick)

En bref

La diffusion d’une espèce dans un milieu : de l’humidité dans un béton, de l’hydrogène dans un acier, des chlorures dans un enrobage. La concentration est l’inconnue, et l’espèce va des zones riches vers les zones pauvres.

species nomme l’espèce : toutes ses grandeurs portent ce nom (c_H2, j_H2, D_H2), si bien que deux espèces vivent sur le même maillage sans se confondre. Ce qui appartient au milieu ne porte pas le nom de l’espèce : le coefficient de stockage poro et les axes matériau. symmetry choisit la forme de la diffusivité, comme pour la conduction.

L’équation est celle de la conduction, mais la physique est distincte : son tag est diffusion, et non thermal, pour qu’un problème qui couple les deux puisse prendre l’une sans l’autre.

Théorie

La première loi de Fick relie le flux de l’espèce au gradient de concentration, et la conservation de l’espèce donne, avec un coefficient de stockage \(\varphi\) (la porosité pour une espèce qui diffuse dans un solide poreux) :

\[ \mathbf j = -\mathsf D\,\nabla c , \qquad \varphi\,\frac{\partial c}{\partial t} - \nabla\cdot(\mathsf D\,\nabla c) = 0 . \]

C’est l’équation de la chaleur avec \(c\) pour \(T\), \(\mathsf D\) pour \(\mathbf K\) et \(\varphi\) pour \(\rho\,c_p\). Sa forme faible et sa discrétisation sont les mêmes :

\[ K_{ij} = \int_e \nabla N_i^{\mathsf T}\,\mathsf D\,\nabla N_j\,d\Omega , \qquad C_{ij} = \int_e \varphi\,N_i\,N_j\,d\Omega . \]

La diffusivité \(\mathsf D\) est un scalaire (isotropic) ou un tenseur orienté par les axes matériau (orthotropic, anisotropic), construit comme la conductivité.

Variables

Pour l’espèce H2 :

La sortie est 𝖣·∇c, le flux de la forme faible ; le flux de Fick est son opposé. Elle est nommée d’après la duale de l’espèce, si bien qu’un modèle qui porte une conduction et plusieurs espèces garde des flux distincts.

Ce qu’elle produit

Exemple

Une barre de diffusivité 0,5 et de longueur 2, où un flux d’hydrogène de 1 entre à gauche, la concentration étant nulle à droite. Le multiplicateur de l’appui est le flux qui sort, égal au flux qui entre :

import pyrucast

# Une barre [0, 2] : un flux d'hydrogène de 1 entre à gauche, c = 0 à droite.
c = pyrucast.Coords(1)
gauche = c.add_node([0.0])
droite = c.add_node([2.0])
barre = pyrucast.mesh.line(gauche, droite, 4)
fes = pyrucast.FiniteElementSpace(barre)

fick = pyrucast.model.fick(fes, "H2")
print(fick[0].primal_vars(), fick[0].dual_vars())  # ['c_H2'] ['j_H2']
appui = pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([droite]), fick, "c_H2")
modele = fick | appui
materiaux = pyrucast.element_field.material_field(modele, [("D_H2", 0.5)])

source = pyrucast.NodeField(barre, ["j_H2"])
source.unit().set_value(gauche, "j_H2", 1.0)
c_h2 = pyrucast.matrix.stiffness(modele, materiaux).solve(
    source | appui.constraint_rhs([(droite, 0.0)])
)
print(round(c_h2.value(gauche, "c_H2"), 9))  # 4.0 = j·L/D
print(
    round(c_h2.value(appui.multiplier_mesh().node(0, 0, 0), "lambda"), 9)
)  # 1.0 : le flux qui sort

gradient = pyrucast.element_field.strain(modele, c_h2, materiaux)
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.components(), flux[0].value(0, 0, "j_H2_x"))  # ['j_H2_x'] -1.0 = D·∂c/∂x

Deux espèces et une conduction sur le même maillage forment un système à trois blocs indépendants :

melange = (
    pyrucast.model.fick(fes, "H2")
    | pyrucast.model.fick(fes, "O2")
    | pyrucast.model.heat_conduction(fes)
)
print(melange.primal_vars())  # ['c_H2', 'c_O2', 'T']
print(len(melange.filter("diffusion")), len(melange.filter("thermal")))  # 2 1
materiaux = pyrucast.element_field.material_field(
    melange, [("D_H2", 0.5), ("D_O2", 0.1), ("k", 2.0)]
)
K = pyrucast.matrix.stiffness(melange, materiaux)  # trois blocs, sans couplage

Une diffusivité qui dépend de la température se calcule par une évolution appliquée au champ de température, avant chaque assemblage. Un transfert d’espèce à travers une interface entre deux corps est un transfert d’interface.

Côté Rust

Fick (src/models/fick.rs) garde son sous-espace, sa MaterialSymmetry et le nom de l’espèce. Ses noyaux, element_stiffness et element_storage, ont la forme de ceux de la conduction et partagent avec eux la construction du tenseur orienté (symmetry::transport_tensor_from). Les noms se construisent depuis l’espèce : fick::primal_var("H2") rend c_H2, fick::dual_var rend j_H2.

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

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

#[test]
fn fick_line_recovers_the_linear_profile() -> Result<()> {
    const D: f64 = 2.0; // diffusivity
    const J: f64 = 10.0; // injected species flux at x = 0
    const C_IMPOSED: f64 = 1.0; // concentration imposed at 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 : diffusion + Dirichlet c = 1 en x = 1 ──────────────────────
    let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(
        nodes.last().unwrap(),
    ))?);
    let diffusion = model::fick(&fes, SPECIES, MaterialSymmetry::Isotropic)?;
    let dirichlet = model::dirichlet(
        &imposed,
        &diffusion,
        &primal_var(SPECIES),
        Default::default(),
    )?;
    let mult = dirichlet.multiplier_mesh()?.node(0, 0, 0)?.id();
    let model = diffusion.union(&dirichlet)?;

    // ── Matériau : diffusivité uniforme ────────────────────────────────────
    let materials =
        pyrucast::ops::element_field::material_field(&model, &[(&format!("D_{SPECIES}"), D)])?;

    // ── Loading: flux J at x = 0, imposed concentration at the multiplier
    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.to_string(), dual_var(SPECIES)])?;
    rhs.set_value(node0, &dual_var(SPECIES), J)?;
    rhs.set_value(mult, IMPOSED, C_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 profile c(x) = 1 + (J/D)(1 − x) ───────
    let tol = 1e-10;
    for (i, node) in nodes.iter().enumerate() {
        let x = i as f64 * h;
        let expected = C_IMPOSED + (J / D) * (1.0 - x);
        let got = solution.value(node.id(), &primal_var(SPECIES))?;
        assert!(
            (got - expected).abs() < tol,
            "c(x={x}) : {got} ≠ {expected}"
        );
    }
    // Mass balance: the reaction at the imposed edge balances the injected flux.
    let reaction = solution.value(mult, MULTIPLIER)?;
    assert!((reaction - J).abs() < tol, "réaction λ : {reaction} ≠ {J}");
    Ok(())
}