Élasticité linéaire

En bref

Un solide qui se déforme réversiblement, proportionnellement à la charge (loi de Hooke), en petites déformations.

kinematics est l’hypothèse cinématique (voir Mécanique), symmetry la symétrie du matériau (voir Orthotropie), strain_measure la mesure de déformation (voir Grandes transformations). Cette page décrit le cas isotrope en petites déformations.

Un maillage de bord (des segments en 2-D) est refusé comme domaine : \(B\) y serait bâti sur le seul gradient tangent, et la raideur n’aurait aucune rigidité dans la direction normale.

Théorie

Sur le domaine \(\Omega\), avec \(b\) les forces de volume :

\[ \nabla\cdot\sigma + b = 0, \qquad \sigma = \mathbb D : \varepsilon, \qquad \varepsilon = \tfrac12\big(\nabla u + \nabla u^{\mathsf T}\big). \]

La forme faible : pour tout déplacement virtuel \(v\),

\[ \int_\Omega \varepsilon(v) : \mathbb D : \varepsilon(u)\,d\Omega = \int_\Omega v\cdot b\,d\Omega + \int_{\Gamma_N} v\cdot t\,d\Gamma , \]

où \(t\) est la traction imposée sur le bord \(\Gamma_N\).

Avec \(u \approx \sum_i N_i\,u_i\), la déformation s’écrit \(\varepsilon = B\,u_e\). En 2-D, \(\varepsilon = [\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy}]\) et le bloc du nœud \(i\) est

\[ B_i = \begin{bmatrix} \partial_x N_i & 0 \\ 0 & \partial_y N_i \\ \partial_y N_i & \partial_x N_i \end{bmatrix}; \]

en 3-D, \(\varepsilon = [\varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz}, \gamma_{yz}, \gamma_{xz}, \gamma_{xy}]\) avec un bloc de six lignes de même forme. La raideur et la masse de l’élément sont

\[ K_e = \int_{\Omega_e} B^{\mathsf T} D\,B\,d\Omega , \qquad M_e = \int_{\Omega_e} \rho\,N^{\mathsf T} N\,d\Omega , \]

intégrées aux points de Gauss, où le matériau est lu.

La matrice de Hooke \(D\), avec \(E\) le module d’Young et \(\nu\) le coefficient de Poisson :

• plane_stress (\(\sigma_{zz} = 0\)), avec \(c = E/(1-\nu^2)\) : \(D = c\,[[1, \nu, 0], [\nu, 1, 0], [0, 0, (1-\nu)/2]]\) ;
• plane_strain (\(\varepsilon_{zz} = 0\)), avec \(c = E/((1+\nu)(1-2\nu))\) : \(D = c\,[[1-\nu, \nu, 0], [\nu, 1-\nu, 0], [0, 0, (1-2\nu)/2]]\) ;
• full_3d, même \(c\) : le bloc normal \(3\times3\) de \(c(1-\nu)\) sur la diagonale et \(c\nu\) ailleurs, et \(G = c\,(1-2\nu)/2\) sur les trois cisaillements ;
• axisymmetric, même \(c\) : le bloc normal isotrope sur \([rr, zz, \theta\theta]\) et \(G\) sur \(rz\).

Axisymétrie. Un solide de révolution se maille dans son plan méridien \((r, z)\). Deux choses changent. La mesure d’intégration devient \(2\pi r\,|J|\), pour toutes les intégrales (raideur, masse, forces) ; c’est la géométrie qui l’apporte, pas la physique. Et la déformation orthoradiale \(\varepsilon_{\theta\theta} = u_r / r\) ajoute une ligne à \(B\), \([N_i/r,\; 0]\) ; elle est rangée dans eps_zz. Les points de Gauss étant intérieurs aux mailles, \(r > 0\) même pour une maille qui touche l’axe.

Dilatation thermique. Une température impose une déformation libre \(\varepsilon_{th} = \alpha\,(T - T_\text{ref})\), sans effet de la mécanique sur la thermique. La raideur ne change pas ; le terme thermique agit sur la contrainte et le second membre :

\[ \sigma = D : (\varepsilon(u) - \varepsilon_{th}), \qquad f_{th} = \int_\Omega B^{\mathsf T} D\,\varepsilon_{th}\,d\Omega . \]

alpha est une composante facultative du matériau ; pyrucast.element_field.thermal_strain(temperature, materiaux, fespace, t_ref) calcule \(\varepsilon_{th}\) aux points de Gauss depuis une température portée aux points de Gauss par interp_to_gauss ; integrate_behavior puis internal_forces en tirent la charge thermique.

Variables

Isotrope, en petites déformations, selon la cinématique :

En axisymétrie, eps_zz et sigma_zz sont les composantes orthoradiales.

Ce qu’elle produit

La raideur géométrique est décrite avec les opérateurs d’assemblage.

Exemple

Une plaque carrée en contraintes planes, tenue en u_x sur son bord gauche et en u_y à son coin inférieur gauche, tirée par une traction \(S = 100\) sur son bord droit. La solution est exacte : \(u_x = S\,x/E\), \(u_y = -\nu\,S\,y/E\) :

gauche, droite, carre = _carre()  # un carré [0, 1]² de 2×2 QUA4
plaque = pyrucast.model.elasticity(pyrucast.FiniteElementSpace(carre), "plane_stress")
coin = gauche.node(0, 0, 0)
modele = (
    plaque
    | pyrucast.model.dirichlet(gauche.to_poi1(), plaque, "u_x")
    | pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([coin]), plaque, "u_y")
    | pyrucast.model.flux(pyrucast.FiniteElementSpace(droite), plaque, "f_x")
)
materiaux = pyrucast.element_field.material_field(
    modele,
    [("E", 200e3), ("nu", 0.3), ("phi_f_x", 100.0)],  # traction S = 100 à droite
)
u = pyrucast.matrix.stiffness(modele, materiaux).solve(
    pyrucast.node_field.external_forces(modele, materiaux)
)
haut_droit = droite.node(0, 1, 1)
print(round(u.value(haut_droit, "u_x"), 12))  # 0.0005 = S/E
print(round(u.value(haut_droit, "u_y"), 12))  # -0.00015 = −ν·S/E

eps = pyrucast.element_field.strain(modele, u, materiaux)
repos = pyrucast.element_field.initial_state(modele, materiaux)
sigma = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, repos)
print(sigma.components())  # ['sigma_xx', 'sigma_yy', 'sigma_xy']
print(round(sigma[0].value(0, 0, "sigma_xx"), 9))  # 100.0

Côté Rust

Elasticity (src/models/elasticity.rs) délègue sa cinématique à un Continuum (src/models/continuum/), qui porte \(B\) pour chaque hypothèse et le noyau de forces internes ; la matrice \(D\) vient de src/models/symmetry.rs, l’isotropie gardant ses formes fermées.

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::model::Model;
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::symmetry::MaterialSymmetry;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::StrainMeasure;
use pyrucast::ops::model;
use pyrucast::ops::solver::solve;
use pyrucast::Result;

#[test]
fn elasticity_unit_square_uniaxial_tension() -> Result<()> {
    const E: f64 = 210.0; // Young's modulus
    const NU: f64 = 0.3; // Poisson's ratio
    const S: f64 = 2.0; // traction on the right edge
    const N: usize = 2; // N×N QUA4 grid
    let h = 1.0 / N as f64;

    // ── QUA4 mesh on [0,1]² ────────────────────────────────────────────────
    let coords = Handle::new(Coords::new(2)?);
    let idx = |i: usize, j: usize| j * (N + 1) + i;
    let mut grid: Vec<Node> = Vec::new();
    for j in 0..=N {
        for i in 0..=N {
            grid.push(Node::create_in(
                coords.clone(),
                &[i as f64 * h, j as f64 * h],
            )?);
        }
    }
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
    for j in 0..N {
        for i in 0..N {
            mesh.add_cell(&[
                grid[idx(i, j)].id(),
                grid[idx(i + 1, j)].id(),
                grid[idx(i + 1, j + 1)].id(),
                grid[idx(i, j + 1)].id(),
            ])?;
        }
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : élasticité plane stress + appuis (rollers) ────────────────
    // The support receives the model it constrains: the dual is read there, and
    // contraints s'y vérifient.
    let roller = |target: &Model, nodes: &[Node], var: &str| -> Result<Model> {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(nodes)?);
        model::dirichlet(&imposed, target, var, Default::default())
    };
    let left: Vec<Node> = (0..=N).map(|j| grid[idx(0, j)].clone()).collect();
    let bottom: Vec<Node> = (0..=N).map(|i| grid[idx(i, 0)].clone()).collect();
    let mut model = model::elasticity(
        &fes,
        Kinematics::PlaneStress,
        MaterialSymmetry::Isotropic,
        StrainMeasure::Linear,
    )?;
    model = model.union(&roller(&model, &left, "u_x")?)?;
    model = model.union(&roller(&model, &bottom, "u_y")?)?;

    // ── Loading: traction S on the right edge (consistent nodal loads, on the
    //    f_x component). The load is a sub-model: it joins the model, its density
    //    the material. ──────────────────────────────────────────────────────
    let mut right_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for j in 0..N {
        right_edge.add_cell(&[grid[idx(N, j)].id(), grid[idx(N, j + 1)].id()])?;
    }
    let right_fes = FiniteElementSpace::lagrange1(&right_edge)?;
    model = model.union(&model::flux(&right_fes, &model, "f_x")?)?;

    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[("E", E), ("nu", NU), ("phi_f_x", S)],
    )?;
    let rhs = pyrucast::ops::node_field::external_forces(&model, &materials)?;

    // ── Assemblage + résolution ────────────────────────────────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical u_x = (S/E)·x, u_y = −(ν S/E)·y ───────
    let tol = 1e-10;
    for j in 0..=N {
        for i in 0..=N {
            let (x, y) = (i as f64 * h, j as f64 * h);
            let ux = solution.value(grid[idx(i, j)].id(), "u_x")?;
            let uy = solution.value(grid[idx(i, j)].id(), "u_y")?;
            assert!((ux - S / E * x).abs() < tol, "u_x({x},{y}) = {ux}");
            assert!((uy + NU * S / E * y).abs() < tol, "u_y({x},{y}) = {uy}");
        }
    }
    Ok(())
}

Sur un solide de révolution, tests/axisymmetric.rs vérifie la solution de Lamé d’un cylindre épais : avec QUA8, l’erreur sur les déplacements nodaux décroît à l’ordre 4 environ et celle des contraintes à l’ordre 2, sans jamais s’annuler, \(1/r\) n’étant pas un polynôme.

La barre chauffée vérifie la dilatation thermique : encastrée aux deux bouts, \(\sigma_{xx} = -E\,\alpha\,\Delta T\) ; libre, une dilatation sans contrainte :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::field::Field;
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::models::tensor::Kinematics;
use pyrucast::ops::element_field::{deformation, interp_to_gauss, thermal_strain, StrainMeasure};
use pyrucast::ops::model;
use pyrucast::ops::node_field::internal_forces;
use pyrucast::ops::solver::solve;
use pyrucast::Result;

#[test]
fn thermoelastic_constrained_bar_stress() -> Result<()> {
    const E: f64 = 210_000.0;
    const NU: f64 = 0.3;
    const ALPHA: f64 = 1e-5;
    const T_REF: f64 = 20.0;
    const DT: f64 = 100.0;
    const NX: usize = 4;
    const NY: usize = 2;
    const L: f64 = 4.0;
    const H: f64 = 1.0;
    let (hx, hy) = (L / NX as f64, H / NY as f64);

    // ── QUA4 mesh on [0,L]×[0,H] ───────────────────────────────────────────
    let coords = Handle::new(Coords::new(2)?);
    let idx = |i: usize, j: usize| j * (NX + 1) + i;
    let mut grid: Vec<Node> = Vec::new();
    for j in 0..=NY {
        for i in 0..=NX {
            grid.push(Node::create_in(
                coords.clone(),
                &[i as f64 * hx, j as f64 * hy],
            )?);
        }
    }
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
    for j in 0..NY {
        for i in 0..NX {
            mesh.add_cell(&[
                grid[idx(i, j)].id(),
                grid[idx(i + 1, j)].id(),
                grid[idx(i + 1, j + 1)].id(),
                grid[idx(i, j + 1)].id(),
            ])?;
        }
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : élasticité + deux bords en x encastrés + appui u_y en bas ──
    let clamp = |target: &Model, nodes: &[Node], var: &str| -> Result<Model> {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(nodes)?);
        model::dirichlet(&imposed, target, var, Default::default())
    };
    let left: Vec<Node> = (0..=NY).map(|j| grid[idx(0, j)].clone()).collect();
    let right: Vec<Node> = (0..=NY).map(|j| grid[idx(NX, j)].clone()).collect();
    let bottom: Vec<Node> = (0..=NX).map(|i| grid[idx(i, 0)].clone()).collect();
    let mut model = model::elasticity(
        &fes,
        Kinematics::PlaneStress,
        MaterialSymmetry::Isotropic,
        StrainMeasure::Linear,
    )?;
    model = model.union(&clamp(&model, &left, "u_x")?)?;
    model = model.union(&clamp(&model, &right, "u_x")?)?;
    model = model.union(&clamp(&model, &bottom, "u_y")?)?;

    // `alpha` supplied through the material field — an optional elastic component.
    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[("E", E), ("nu", NU), ("alpha", ALPHA)],
    )?;

    // ── Imposed temperature T = T_ref + ΔT everywhere, carried at the Gauss points
    let support = Handle::new(SubMesh::poi1_from_nodes(&grid)?);
    let mut t_nodal = SubNodeField::from_poi1(&support, vec!["T".into()])?;
    for n in &grid {
        t_nodal.set_value(n.id(), "T", T_REF + DT)?;
    }
    let t_elem = interp_to_gauss(&NodeField::from_sub(t_nodal), &fes)?;

    // ── Charge thermique f_th = ∫ Bᵀ D ε_th (BSIG de σ_th = D:ε_th) ─────────
    let eps_th = thermal_strain(&t_elem, &materials, &fes, T_REF, StrainMeasure::Linear)?;
    let sig_th = pyrucast::ops::element_field::behavior::integrate(
        &model,
        &eps_th,
        &pyrucast::ops::element_field::behavior::initial_state(&model, &materials)?,
        &materials,
        None,
    )?;
    let f_th = internal_forces(&model, &sig_th, &NodeField::empty(), &materials)?;

    // ── Assemblage + résolution ────────────────────────────────────────────
    let solution = solve(
        &pyrucast::ops::matrix::stiffness(&model, &materials)?,
        &f_th,
    )?;

    // ── Déplacement propre (u_x, u_y) puis σ = D:(ε(u) − ε_th) ─────────────
    let disp_support = Handle::new(SubMesh::poi1_from_nodes(&grid)?);
    let mut disp = SubNodeField::from_poi1(&disp_support, vec!["u_x".into(), "u_y".into()])?;
    for n in &grid {
        disp.set_value(n.id(), "u_x", solution.value(n.id(), "u_x")?)?;
        disp.set_value(n.id(), "u_y", solution.value(n.id(), "u_y")?)?;
    }
    let eps = deformation(&NodeField::from_sub(disp), &fes, StrainMeasure::Linear)?;
    let eps_mech = eps.merge_field(&eps_th, |a, b| a - b)?;
    let rest = pyrucast::ops::element_field::behavior::initial_state(&model, &materials)?;
    let sigma = pyrucast::ops::element_field::behavior::integrate(
        &model, &eps_mech, &rest, &materials, None,
    )?;

    // ── Vérification : σ_xx = −E·α·ΔT, σ_yy = 0 ────────────────────────────
    let expected = -E * ALPHA * DT;
    let tol = 1e-6 * expected.abs();
    let sub = sigma.get(0)?.read();
    for cell in 0..sub.cell_count() {
        for g in 0..sub.gauss_count() {
            assert!((sub.value(cell, g, "sigma_xx")? - expected).abs() < tol);
            assert!(sub.value(cell, g, "sigma_yy")?.abs() < tol);
        }
    }
    Ok(())
}