Diffusion (loi de Fick)

Introduction

La diffusion d’une espèce dans un milieu — humidité dans un béton, hydrogène dans un acier, chlorures dans un enrobage — obéit à la loi de Fick, dont l’opérateur est celui de la conduction thermique.

C’est pourtant une physique distincte, et pyrucast la traite comme telle. La variable primale est la concentration c, la duale le flux de matière j, et sa nature déclarée est Physics::Diffusion. Partager un opérateur n’est pas partager une physique : dans un problème couplé thermo-diffusif, on doit pouvoir écrire model.filter("diffusion") et n’obtenir que la partie diffusive, sans traîner la thermique avec.

Le modèle vit sur n’importe quel espace EF volumique (2-D ou 3-D, linéaire ou quadratique), avec un degré de liberté scalaire par nœud.

Équations continues résolues

La première loi de Fick relie le flux au gradient de concentration :

\[ \mathbf j = -\,\mathsf D\,\nabla c \]

et la conservation de l’espèce, en régime transitoire avec un coefficient de stockage \( \varphi \) (la porosité, pour une espèce diffusant dans un solide poreux) :

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

En stationnaire c’est l’équation de Laplace pondérée par \( \mathsf D \). C’est la même équation que la conduction thermique, à un changement de noms près (\( c \leftrightarrow T \), \( \mathsf D \leftrightarrow k \), \( \varphi \leftrightarrow \rho c_p \)) — d’où un modèle qui partage la totalité du noyau de conduction thermique, et n’en diffère que par ses variables et sa nature physique.

La forme faible, après intégration par parties, s’écrit : trouver c tel que pour tout δc admissible,

\[ \int_\Omega \varphi\,\delta c\,\dot c\; d\Omega

• \int_\Omega \nabla \delta c \cdot \mathsf D\,\nabla c\; d\Omega = -\int_{\partial\Omega} \delta c\;\mathbf j\!\cdot\!\mathbf n\; d\Gamma . \]

Forme discrétisée

\[ K_{ij} = \int_\Omega \nabla N_i^\top\,\mathsf D\;\nabla N_j\; d\Omega \quad \text{(rigidité de diffusion — Cast3M \texttt{COND})}, \] \[ C_{ij} = \int_\Omega \varphi\,N_i\,N_j\; d\Omega \quad \text{(stockage — Cast3M \texttt{CAPA})}. \]

D est un tenseur, dont le cas isotrope D = D·I redonne le produit scalaire habituel ∇N_i · ∇N_j. Les trois symétries matériau décrites au chapitre Élasticité orthotrope s’appliquent identiquement, avec un tenseur d’ordre 2 au lieu de 4 :

Le repère est donné par les vecteurs V1X, V1Y (2-D) ou V1X…V1Z, V2X…V2Z (3-D), exactement comme en mécanique.

Variables et matériau

Le comportement (COMP) rend le flux sous forme faible D·∇c, en composantes j_x, j_y(, j_z). Comme en thermique, c’est l’opposé du flux physique de Fick : ce choix garantit ∫ Bᵀ·j = K·c, donc l’accord entre le comportement et la rigidité dans le cas linéaire. Les composantes sont nommées d’après la variable duale (j_*) et non flux_*, afin qu’un modèle portant à la fois conduction et diffusion garde deux champs de flux non ambigus.

L’entrée du comportement est le gradient grad_c_x, …, tel que le produit l’opérateur gradient sur un champ dont la composante est c.

Mise en donnée (Rust, testé)

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::Physics;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::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 multiplier = mesh::barycenter(&imposed)?;
    let mult = multiplier.node(0, 0, 0)?.id();

    let diffusion = model::fick(&fes, SPECIES)?;
    let dirichlet = model::dirichlet(
        &diffusion,
        &primal_var(SPECIES),
        &imposed,
        &multiplier,
        Default::default(),
    )?;
    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![
            format!("imposed_{}", primal_var(SPECIES)),
            dual_var(SPECIES),
        ],
    )?;
    rhs.set_value(node0, &dual_var(SPECIES), J)?;
    rhs.set_value(mult, &format!("imposed_{}", primal_var(SPECIES)), 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, &format!("lambda_{}", primal_var(SPECIES)))?;
    assert!((reaction - J).abs() < tol, "réaction λ : {reaction} ≠ {J}");
    Ok(())
}

Exemple Python

"""1-D Fick diffusion — a bar fed with a species, compared with the analytical solution.

Problème
--------
On the segment [0, 1]:

  * en x = 0 : un **flux d'espèce** imposé ``J`` (Neumann) ;
  * at x = 1: an **imposed concentration** ``c = 1`` (Dirichlet).

In the steady regime without a volume source the profile is linear ::

    c(x) = 1 + (J / D) * (1 - x)

and the Lagrange multiplier at the imposed node is exactly ``J``: everything
entering at x = 0 leaves at x = 1 (mass balance).

The operator is the thermal conduction one; what changes is the **physics**.
The primal is the concentration ``c``, the dual the flux ``j``, and the kind
declared is ``"diffusion"`` — so that a coupled thermo-diffusive model splits
with ``model.filter(...)``, which the second part
de l'exemple montre.

This is the Python equivalent of the Rust integration test ``tests/fick.rs``.

Lancement
---------
Once the extension is built in the venv ::

    maturin develop --features extension-module
    python examples/diffusion_1d.py
"""

import pyrucast

# ── Problem data ────────────────────────────────────────────────────────────
SPECIES = "H2"  # the diffusing species — every name carries it
D = 2.0  # diffusivité
J = 10.0  # flux d'espèce injecté en x = 0
C_IMPOSED = 1.0  # concentration imposée en x = 1
N_ELEMS = 4
K = 5.0  # thermal conductivity, for the coupled part


def ligne(n_elems):
    """A line of ``n_elems`` SEG2 on [0, 1], with its nodes."""
    c = pyrucast.Coords(1)
    h = 1.0 / n_elems
    nodes = [c.add_node([i * h]) for i in range(n_elems + 1)]
    mesh = pyrucast.Mesh(c, "SEG2")
    for i in range(n_elems):
        mesh.unit().add_cell([nodes[i], nodes[i + 1]])
    return c, nodes, pyrucast.FiniteElementSpace(mesh), h


def profil_stationnaire() -> None:
    c, nodes, fes, h = ligne(N_ELEMS)

    # ── Modèle : diffusion + concentration imposée en x = 1 ──────────────────
    imposed = pyrucast.Mesh(c, "POI1")
    imposed.unit().add_cell([nodes[-1]])
    multiplier = pyrucast.mesh.barycenter(imposed)
    mult = multiplier.node(0, 0, 0)

    cible = pyrucast.model.fick(fes, SPECIES)

    model = cible | pyrucast.model.dirichlet(cible, f"c_{SPECIES}", imposed, multiplier)
    materials = pyrucast.element_field.material_field(model, [(f"D_{SPECIES}", D)])

    # ── Loading: flux J at x = 0, imposed value at the multiplier ────────────
    load = pyrucast.Mesh(c, "POI1")
    load.unit().add_cell([nodes[0]])
    load.unit().add_cell([mult])
    rhs = pyrucast.NodeField(load, [f"imposed_c_{SPECIES}", f"j_{SPECIES}"])
    rhs[0].set_value(nodes[0], f"j_{SPECIES}", J)
    rhs[0].set_value(mult, f"imposed_c_{SPECIES}", C_IMPOSED)

    # ── Assemblage + résolution ─────────────────────────────────────────────
    stiffness = pyrucast.matrix.stiffness(model, materials)
    solution = pyrucast.solver.solve(stiffness, rhs)

    print("Diffusion de Fick 1-D")
    print(f"  D = {D}, flux injecté J = {J}, c(1) = {C_IMPOSED}")
    print()
    print("     x      c calculé    c analytique")
    print("  " + "-" * 36)
    for i, node in enumerate(nodes):
        x = i * h
        attendu = C_IMPOSED + (J / D) * (1.0 - x)
        obtenu = solution.value(node, f"c_{SPECIES}")
        print(f"  {x:5.3f}   {obtenu:10.6f}   {attendu:12.6f}")
        assert abs(obtenu - attendu) < 1e-10

    reaction = solution.value(mult, f"lambda_c_{SPECIES}")
    print()
    print(f"  Bilan de matière : réaction = {reaction:.6f}, flux injecté = {J}")
    assert abs(reaction - J) < 1e-10


def couplage_avec_la_thermique() -> None:
    """Diffusion and conduction on the same mesh: two distinct physics."""
    _c, _nodes, fes, _h = ligne(3)
    model = pyrucast.model.fick(fes, SPECIES) | pyrucast.model.heat_conduction(fes)

    # A single material field carries both sets: the assembler resolves each zone
    # through the components its physics requires (`D` here, `k` there).
    materials = pyrucast.element_field.material_field(
        model, [(f"D_{SPECIES}", D), ("k", K)]
    )
    pyrucast.matrix.stiffness(model, materials)

    print()
    print("Modèle couplé diffusion + thermique")
    print(f"  sous-modèles          : {len(model)}")
    print(f"  filter('diffusion')   : {len(model.filter('diffusion'))}")
    print(f"  filter('thermal')     : {len(model.filter('thermal'))}")
    print(f"  filter('mechanical')  : {len(model.filter('mechanical'))}")
    assert len(model.filter("diffusion")) == 1
    assert len(model.filter("thermal")) == 1
    assert len(model.filter("mechanical")) == 0


def main() -> None:
    profil_stationnaire()
    couplage_avec_la_thermique()


if __name__ == "__main__":
    main()

Compléments

Coexister avec la thermique

Les deux physiques peuvent vivre sur le même maillage sans se gêner :

model = pyrucast.model.fick(fes, "H2") | pyrucast.model.heat_conduction(fes)
materials = pyrucast.element_field.material_field(model, [("D_H2", 2.0), ("k", 5.0)])
k = pyrucast.matrix.stiffness(model, materials)

len(model.filter("diffusion"))  # 1
len(model.filter("thermal"))  # 1

Un seul champ matériau porte les deux jeux de coefficients. L’assembleur résout la zone de chaque physique par les composantes qu’elle exige (D ici, k là) — il n’y a rien à consolider à la main. Et parce que les deux natures sont distinctes, filter les sépare de nouveau après coup, aussi bien sur le modèle que sur la matrice assemblée.

Les degrés de liberté restent séparés (c d’un côté, T de l’autre) : le système est bloc-diagonal. Un vrai couplage — une diffusivité fonction de la température, ou une thermodiffusion — se pilote depuis Python, en réassemblant la partie diffusive à chaque pas avec un champ matériau recalculé.

Transfert à travers une interface

Deux corps qui se touchent ne partagent pas forcément leurs nœuds. Un contact imparfait, un revêtement, un joint, une membrane laissent le champ sauter à la traversée, tandis qu’un flux la franchit proportionnellement à ce saut :

\[ j\cdot n = h\,\big(c_1 - c_2\big) \]

h_c_H2 est le coefficient de transfert (son inverse est la résistance de contact) : un par grandeur transférée, nommé d’après elle.

Ce modèle n’a rien de diffusif non plus : on lui passe [("T", "q")] et une conduction pour cible pour une résistance de contact, les couples de déplacement et une élasticité pour un joint collé de raideur finie — la nature vient de la cible. La loi commune, sa structure en quatre blocs dont deux hors-diagonale, l’exigence de conformité des deux côtés et le critère qui départage un échange d’une contrainte MPC sont dans Échanges.

corps = pyrucast.model.fick(gauche, "H2") | pyrucast.model.fick(droite, "H2")
model = corps | pyrucast.model.interface_transfer(
    face_gauche, face_droite, corps, [("c_H2", "j_H2")]
)
materials = pyrucast.element_field.material_field(
    model, [("D_H2", 2.0), ("h_c_H2", 5.0)]
)

Ce que ça vaut comme vérification

Deux carrés côte à côte, un flux q injecté d’un côté, la concentration imposée de l’autre : le profil est linéaire par morceaux avec une chute q/D dans chaque carré et un saut q/h à l’interface. C’est ce saut qui distingue une interface d’un nœud partagé, et il est porté entièrement par les blocs hors-diagonale. Quand h → ∞, le saut s’efface et l’on retrouve le corps continu.

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

Régime transitoire

La matrice de stockage s’assemble avec matrix.mass(...), qui exige alors la composante poro. L’intégration en temps est orchestrée en Python, comme pour la thermique transitoire — le noyau Rust fournit K et C, pas la boucle.

Bilan de matière

Comme en thermique, le multiplicateur de Lagrange d’une concentration imposée est le flux d’espèce qui traverse la frontière. C’est la vérification la plus directe d’un calcul de diffusion, et c’est ce que contrôle le test d’intégration ci-dessus : la réaction au bord imposé égale exactement le flux injecté à l’autre bout.