Poutre de Timoshenko

Poutre déformable en cisaillement : élément SEG2, dans la configuration que lui donne la dimension du maillage. C’est une seule physique — elle remplace les anciens Model.frame (portique plan) et Model.frame3d (cadre spatial), qui en étaient les cas 2-D et 3-D.

Il n’y a rien à choisir : tout ce qui distingue les trois — nombre et noms des DDL, jeu matériau, efforts rendus, présence d’un terme axial, d’une torsion, d’une rotation vers les axes globaux — découle de la dimension. On relit les noms obtenus par model.primal_vars().

Équations continues résolues

La section reste plane mais non normale à l’axe déformé : la rotation θ est un champ indépendant, et la distorsion γ = w' − θ une déformation à part entière. C’est toute la différence avec Euler-Bernoulli, où θ = w' et où γ n’existe pas.

• cinématique : courbure \( \kappa = \theta’ \), distorsion \( \gamma = w’ - \theta \) ;
• efforts : \( M = EI,\theta’ \), \( V = G A_s (w’ - \theta) \) ;
• équilibre : \( V’ + q = 0 \), \( M’ - V = 0 \).

Ces deux équations sont du second ordre, là où Bernoulli en a une seule du quatrième. C’est la contrepartie de l’hypothèse cinématique : libérer θ de w' abaisse l’ordre de l’équation, et abaisse avec lui l’exigence de continuité — le C⁰ suffit là où Bernoulli réclame du C¹.

Forme discrétisée — l’élément exact

L’élément assemblé est la solution exacte de ces deux équations sur une travée libre d’efforts répartis. Ses fonctions de forme sont cubiques en w et quadratiques en θ, et elles portent le matériau par

\[ \Phi = \frac{12,E I}{G A_s L^2}, \]

le rapport des souplesses de flexion et de cisaillement. La flexion s’écrit alors en forme fermée :

\[ K_b = \frac{EI}{L^3(1+\Phi)} \begin{bmatrix} 12 & 6L & -12 & 6L \\ 6L & (4+\Phi)L^2 & -6L & (2-\Phi)L^2 \\ -12 & -6L & 12 & -6L \\ 6L & (2-\Phi)L^2 & -6L & (4+\Phi)L^2 \end{bmatrix}. \]

L’élément est exact aux nœuds pour des charges d’extrémité : un élément par barre suffit. On lit directement sur cette matrice que la raideur en flèche est la combinaison en série des deux souplesses,

\[ K_{ww} = \frac{12EI}{L^3(1+\Phi)} = \frac{1}{\dfrac{L^3}{12EI} + \dfrac{L}{G A_s}}, \]

— on fléchit le tronçon et on le cisaille, les deux cèdent l’un après l’autre. Et \( \Phi = 0 \) redonne terme pour terme la matrice d’Euler-Bernoulli : « Bernoulli est la limite sans cisaillement » est une propriété vérifiée par un test, pas une phrase.

L’espace EF ne porte aucune base

Ces fonctions de forme dépendent du matériau par \( \Phi \). Aucun espace éléments finis ne peut donc les tabuler — il tabule par type d’élément, pas par maille. L’espace déclare en conséquence MODEL_EMBEDDED : la formulation possède son interpolation, et le dit.

fes = pyrucast.FiniteElementSpace(maillage, interpolation="MODEL_EMBEDDED")
poutre = pyrucast.model.timoshenko(fes)

Ce que remplace cet élément. La version précédente était linéaire, à cisaillement sous-intégré : elle convergeait au raffinement au lieu d’être exacte, et déclarait une interpolation de Lagrange qu’elle utilisait réellement. Le portique 2-D était dans ce cas, le cadre 3-D employait déjà la forme exacte — deux modèles frères, deux théories discrètes. Ils n’en font plus qu’une.

Variables et matériau

Voir le tableau d’ouverture. rho est facultatif, exigé par la seule matrice de masse ; en configuration 1-D l’aire pleine A l’est aussi, la rigidité n’utilisant que l’aire de cisaillement.

Le comportement (COMP) rend les efforts de section par une loi linéaire, à partir des déformations généralisées produites par beam_deformation, un opérateur pour les trois configurations.

Mise en donnée (Rust, testé)

Console encastrée, charge transverse P au bout libre ; solution analytique w = P·L³/(3EI) + P·L/(G·A_s) — les deux souplesses, en série.

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Interpolation, Node};
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::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn timoshenko_cantilever_converges_without_locking() -> Result<()> {
    const E: f64 = 1.0;
    const I: f64 = 1.0; // E·I = 1
    const G: f64 = 30.0;
    const A_S: f64 = 1.0; // G·A_s = 30 (slender ⇒ shear locking would be severe)
    const L: f64 = 1.0;
    const P: f64 = 1.0; // transverse tip load
    const N: usize = 40; // beam elements

    // ── Mesh: N SEG2 elements aligned on [0, L] (1-D configuration) ────────
    let coords = Handle::new(Coords::new(1)?);
    let h = L / N as f64;
    let nodes: Vec<Node> = (0..=N)
        .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 {
        mesh.add_cell(&[nodes[i].id(), nodes[i + 1].id()])?;
    }
    let fes = FiniteElementSpace::new(&mesh, Interpolation::ModelEmbedded)?;

    // ── Model: beam + clamping on the left (w = θ = 0) ─────────────────────
    let clamp = |target: &Model, node: &Node, var: &str| -> Result<Model> {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(node))?);
        let multiplier = mesh::barycenter(&imposed)?;
        model::dirichlet(target, var, &imposed, &multiplier, Default::default())
    };
    let mut model = model::timoshenko(&fes)?;
    model = model.union(&clamp(&model, &nodes[0], "w")?)?;
    model = model.union(&clamp(&model, &nodes[0], "theta")?)?;

    // ── Matériau E, I, G, A_s ──────────────────────────────────────────────
    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[("E", E), ("I", I), ("G", G), ("A_s", A_S)],
    )?;

    // ── Loading: transverse force P at the free end (component f_w) ────────
    let mut load_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    load_sm.add_cell(&[nodes[N].id()])?;
    let load_sm = Handle::new(load_sm);
    let mut rhs = SubNodeField::from_poi1(&load_sm, vec!["f_w".into()])?;
    rhs.set_value(nodes[N].id(), "f_w", P)?;
    let rhs = NodeField::from_sub(rhs);

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

    // ── Comparaison : w_tip = P·L³/(3·E·I) + P·L/(G·A_s) ───────────────────
    let w_tip = solution.value(nodes[N].id(), "w")?;
    let analytical = P * L.powi(3) / (3.0 * E * I) + P * L / (G * A_S);
    assert!(
        (w_tip - analytical).abs() < 1e-2 * analytical,
        "w_tip = {w_tip}, analytique {analytical}"
    );
    Ok(())
}

Le portique plan, où l’axial et la flexion se découplent :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Interpolation, Node};
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::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn frame_inclined_cantilever_perpendicular_load() -> Result<()> {
    const E: f64 = 1.0;
    const A: f64 = 1.0;
    const I: f64 = 1.0;
    const G: f64 = 30.0;
    const A_S: f64 = 1.0;
    const L: f64 = 1.0;
    const P: f64 = 1.0;
    const N: usize = 40;

    // Beam direction (45°) and perpendicular.
    let (c, s) = (
        std::f64::consts::FRAC_1_SQRT_2,
        std::f64::consts::FRAC_1_SQRT_2,
    );
    let (px, py) = (-s, c); // unit perpendicular
    let h = L / N as f64;

    // ── Mesh: N SEG2 elements along the 45° direction ──────────────────────
    let coords = Handle::new(Coords::new(2)?);
    let nodes: Vec<Node> = (0..=N)
        .map(|i| Node::create_in(coords.clone(), &[i as f64 * h * c, i as f64 * h * s]))
        .collect::<Result<_>>()?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for i in 0..N {
        mesh.add_cell(&[nodes[i].id(), nodes[i + 1].id()])?;
    }
    let fes = FiniteElementSpace::new(&mesh, Interpolation::ModelEmbedded)?;

    // ── Model: frame + full clamping at the base ───────────────────────────
    let clamp = |target: &Model, node: &Node, var: &str| -> Result<Model> {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(node))?);
        let multiplier = mesh::barycenter(&imposed)?;
        model::dirichlet(target, var, &imposed, &multiplier, Default::default())
    };
    let mut model = model::timoshenko(&fes)?;
    model = model.union(&clamp(&model, &nodes[0], "u_x")?)?;
    model = model.union(&clamp(&model, &nodes[0], "u_y")?)?;
    model = model.union(&clamp(&model, &nodes[0], "r_z")?)?;

    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[("E", E), ("A", A), ("I", I), ("G", G), ("A_s", A_S)],
    )?;

    // ── Loading: force P perpendicular to the beam, at the free end ────────
    let mut load_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    load_sm.add_cell(&[nodes[N].id()])?;
    let load_sm = Handle::new(load_sm);
    let mut rhs = SubNodeField::from_poi1(&load_sm, vec!["f_x".into(), "f_y".into()])?;
    rhs.set_value(nodes[N].id(), "f_x", P * px)?;
    rhs.set_value(nodes[N].id(), "f_y", P * py)?;
    let rhs = NodeField::from_sub(rhs);

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

    // ── Comparison: the tip's displacement = δ·(perpendicular) ─────────────
    let delta = P * L.powi(3) / (3.0 * E * I) + P * L / (G * A_S);
    let ux = solution.value(nodes[N].id(), "u_x")?;
    let uy = solution.value(nodes[N].id(), "u_y")?;
    // Projection onto the perpendicular (= δ) and onto the axis (≈ 0).
    let transverse = ux * px + uy * py;
    let axial = ux * c + uy * s;
    assert!(
        (transverse - delta).abs() < 1e-2 * delta,
        "transverse {transverse} ≠ {delta}"
    );
    assert!(axial.abs() < 1e-6, "déplacement axial {axial} ≈ 0");
    Ok(())
}

Et le cadre spatial, avec sa torsion :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Interpolation, Node};
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::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn frame3d_cantilever_bending_and_torsion() -> Result<()> {
    const E: f64 = 1.0;
    const A: f64 = 1.0;
    const IY: f64 = 1.0;
    const IZ: f64 = 2.0;
    const J: f64 = 1.0;
    const G: f64 = 0.5;
    const ASY: f64 = 10.0;
    const ASZ: f64 = 10.0;
    const L: f64 = 1.0;
    const PY: f64 = 1.0;
    const PZ: f64 = 1.0;
    const MX: f64 = 1.0;
    const N: usize = 2;

    // ── Maillage : N éléments SEG2 le long de l'axe X (config 3-D) ─────────
    let coords = Handle::new(Coords::new(3)?);
    let h = L / N as f64;
    let nodes: Vec<Node> = (0..=N)
        .map(|i| Node::create_in(coords.clone(), &[i as f64 * h, 0.0, 0.0]))
        .collect::<Result<_>>()?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for i in 0..N {
        mesh.add_cell(&[nodes[i].id(), nodes[i + 1].id()])?;
    }
    let fes = FiniteElementSpace::new(&mesh, Interpolation::ModelEmbedded)?;

    // ── Model: 3-D frame + full clamping (6 DOFs) at the base ──────────────
    let clamp = |target: &Model, node: &Node, var: &str| -> Result<Model> {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(node))?);
        let multiplier = mesh::barycenter(&imposed)?;
        model::dirichlet(target, var, &imposed, &multiplier, Default::default())
    };
    let mut model = model::timoshenko(&fes)?;
    for var in ["u_x", "u_y", "u_z", "r_x", "r_y", "r_z"] {
        model = model.union(&clamp(&model, &nodes[0], var)?)?;
    }

    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[
            ("E", E),
            ("A", A),
            ("I_y", IY),
            ("I_z", IZ),
            ("J", J),
            ("G", G),
            ("A_sy", ASY),
            ("A_sz", ASZ),
        ],
    )?;

    // ── Loading: f_y, f_z and m_x at the free end ──────────────────────────
    let mut load_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    load_sm.add_cell(&[nodes[N].id()])?;
    let load_sm = Handle::new(load_sm);
    let mut rhs =
        SubNodeField::from_poi1(&load_sm, vec!["f_y".into(), "f_z".into(), "m_x".into()])?;
    rhs.set_value(nodes[N].id(), "f_y", PY)?;
    rhs.set_value(nodes[N].id(), "f_z", PZ)?;
    rhs.set_value(nodes[N].id(), "m_x", MX)?;
    let rhs = NodeField::from_sub(rhs);

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

    // ── Compared with the analytical solution (exact element ⇒ nodally exact) ─
    let tip = nodes[N].id();
    let uy = PY * L.powi(3) / (3.0 * E * IZ) + PY * L / (G * ASY);
    let uz = PZ * L.powi(3) / (3.0 * E * IY) + PZ * L / (G * ASZ);
    let rx = MX * L / (G * J);
    let tol = 1e-9;
    assert!((solution.value(tip, "u_y")? - uy).abs() < tol, "u_y");
    assert!((solution.value(tip, "u_z")? - uz).abs() < tol, "u_z");
    assert!((solution.value(tip, "r_x")? - rx).abs() < tol, "r_x");
    // Axial DOF stays put (no axial load).
    assert!(solution.value(tip, "u_x")?.abs() < tol, "u_x ≈ 0");
    Ok(())
}

Exemple Python

"""Poutre de Timoshenko — console exacte dès un seul élément.

Physique
--------
Poutre déformable en cisaillement. Cinématique : courbure `κ = θ'`, distorsion
`γ = w' - θ`. Efforts : moment `M = E·I·θ'`, effort tranchant
`V = G·A_s·(w' - θ)`. Équilibre : `dV/dx + q = 0`, `dM/dx - V = 0`.

The assembled element is the **exact solution** of these two equations on a
span free of distributed loads — the closed form parameterized by
`Φ = 12·E·I/(G·A_s·L²)`. Its shape functions therefore depend on the material,
which no finite element space can tabulate: the space declares
`MODEL_EMBEDDED`, that is, the formulation owns its interpolation.

Problème
--------
Clamped cantilever (`w = θ = 0`), transverse load `P` at the free end.
Analytical solution `w = P·L³/(3·E·I) + P·L/(G·A_s)` — both compliances,
cisaillement, **en série**.

The element being exact at the nodes, **one** is enough: refining changes
nothing, which this script checks. (The previous version was linear with
under-integrated shear; it converged towards that value instead of reaching
it, and this example showed its convergence.)

Lancement ::

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

import pyrucast

E, I, G, A_S, L, P = 1.0, 1.0, 30.0, 1.0, 1.0, 1.0


def _clamp(target, node, var):
    imposed = pyrucast.mesh.poi1_from_nodes([node])
    multiplier = pyrucast.mesh.barycenter(imposed)
    return pyrucast.model.dirichlet(target, var, imposed, multiplier)


def tip_deflection(n_elems: int) -> float:
    c = pyrucast.Coords(1)
    base = c.add_node([0.0])
    tip = c.add_node([L])
    mesh = pyrucast.mesh.line(base, tip, n_elems)  # console 1-D (`line`)
    # The basis belongs to the formulation, not to the space: it depends on `Φ`,
    # hence on the material, and is computed cell by cell.
    fes = pyrucast.FiniteElementSpace(mesh, interpolation="MODEL_EMBEDDED")

    model = pyrucast.model.timoshenko(fes)
    model = model | _clamp(model, base, "w")
    model = model | _clamp(model, base, "theta")

    materials = pyrucast.element_field.material_field(
        model, [("E", E), ("I", I), ("G", G), ("A_s", A_S)]
    )

    load = pyrucast.mesh.poi1_from_nodes([tip])
    rhs = pyrucast.NodeField(load, ["f_w"])
    rhs[0].set_value(tip, "f_w", P)

    solution = pyrucast.solver.solve(pyrucast.matrix.stiffness(model, materials), rhs)
    return solution.value(tip, "w")


def main() -> None:
    analytical = P * L**3 / (3.0 * E * I) + P * L / (G * A_S)
    print(f"{'N':>4} {'w_tip':>12} {'err. rel.':>12}")
    for n in (1, 2, 5, 10, 40):
        w = tip_deflection(n)
        print(f"{n:4d} {w:12.6f} {abs(w - analytical) / analytical:12.2e}")
    print(f"\nanalytique  = {analytical:.6f}  (P·L³/3EI + P·L/GA_s)")

    # Exact at the nodes: one element already gives the answer, and refining does
    # not improve it — there is nothing to improve.
    one = tip_deflection(1)
    assert abs(one - analytical) < 1e-12 * analytical, one
    assert abs(tip_deflection(40) - one) < 1e-12 * analytical
    print("OK: exact at the nodes from one element on, refining changes nothing.")


if __name__ == "__main__":
    main()

Compléments

Masse. La masse cohérente est celle du même élément, intégrée des mêmes fonctions de forme que sa rigidité :

\[ M = \int_0^L \rho A\, N_w^\top N_w\, dx

• \int_0^L \rho I\, N_\theta^\top N_\theta\, dx, \]

le second terme étant l’inertie de rotation de la section. Rigidité et masse décrivent enfin une seule poutre. Seuls l’axial et la torsion gardent la forme du champ linéaire (ρL/6)[[2,1],[1,2]], exacte pour ce qu’ils interpolent réellement.

Elle est intégrée, et non recopiée de la table publiée de polynômes en Φ. Cette table est juste, mais un coefficient mal retranscrit sur vingt donnerait une matrice plausible, symétrique et définie positive — décrivant une autre poutre. C’est le mode de défaillance qui avait coûté une tangente fausse plus tôt dans ce projet. Une intégration ne se retranscrit pas : les fonctions de forme sont celles de la rigidité, et quatre points de Gauss rendent la quadrature exacte (l’intégrande est de degré 6).

Ce qui l’épingle : à Φ = 0 elle redonne la table classique ρAL/420 · [156, 22L, 54, −13L ; …], douze nombres que personne ne conteste ; une translation rigide porte exactement ρAL quel que soit Φ ; et le couplage flèche-rotation, absent de la masse linéaire, est bien là.

Rigidité géométrique. Elle demande un effort axial pour raidir la barre : la configuration 1-D, en flexion pure, n’en déclare donc aucune.

Reconstruction des efforts. beam_deformation évalue les déformations à chaque point de Gauss, depuis les fonctions de forme de l’élément — les mêmes que sa rigidité et sa masse. Ce qu’elle rend dit alors la physique :

• la courbure varie linéairement, donc le moment aussi, ce qu’impose M' = V ;
• le cisaillement est constant, ce qu’impose V' = 0 sur une travée non chargée.

L’élément linéaire ne pouvait rendre ni l’un ni l’autre : sa courbure était constante et son cisaillement oscillait, d’où une moyenne pour tout.

L’opérateur exige le matériau, et c’est la signature honnête : Φ en dépend, donc la distribution de courbure aussi. On ne reconstitue pas la courbure d’une poutre sans connaître sa raideur de cisaillement.

Le résidu, et pourquoi Φ est dans l’état. Les forces internes intègrent le transposé du même B (models::beam::b_into), si bien que ∫ Bᵀσ vaut K·u exactement — la loi de section étant linéaire. Mais le noyau qui les calcule reçoit la géométrie et l’état, jamais le matériau : le B d’un continuum est le gradient symétrique et ignore tout module, ce seam n’a donc jamais eu à en porter un. Le comportement rend par conséquent Φ (phi, ou phi_y et phi_z en spatial) à côté des efforts de section. C’est la seule grandeur non conjuguée que ce dépôt garde en état, et elle le mérite : le résidu la relit à chaque itération de Newton.

Le repère local, lui, est déduit automatiquement de la géométrie (référence globale Z, ou Y pour une barre verticale), ce qui convient aux sections symétriques.

Le repère de rotation. Le portique plan nommait sa rotation rz quand tout le reste du dépôt écrivait r_z. La fusion l’a fait sortir immédiatement — une matrice non carrée au solveur — et c’est r_z qui l’emporte.