Poutre de Timoshenko

En bref

Une poutre déformable en cisaillement : la section reste plane mais pas normale à l’axe déformé. Pour une poutre trapue, la flèche due au cisaillement s’ajoute à celle de la flexion.

Comme pour Euler-Bernoulli, la configuration se lit sur la dimension des Coords : flexion pure en 1-D, portique plan en 2-D, cadre spatial en 3-D.

Théorie

La rotation \(\theta\) est un champ indépendant de la pente, et la distorsion \(\gamma = w' - \theta\) une déformation à part entière :

\[ \kappa = \theta', \qquad \gamma = w' - \theta, \qquad M = E I\,\kappa, \qquad V = G A_s\,\gamma, \qquad V' + q = 0, \qquad M' - V = 0 , \]

\(A_s\) étant l’aire de cisaillement. Deux équations du second ordre : une interpolation \(C^0\) suffit.

L’élément est la solution exacte de ces équations sur une travée sans charge répartie : il est exact aux nœuds, et un élément par barre suffit. Ses fonctions de forme, cubiques en \(w\) et quadratiques en \(\theta\), dépendent du matériau par

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

le rapport des souplesses de flexion et de cisaillement. La raideur de flexion est

\[ K_b = \frac{E I}{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}, \]

où la raideur en flèche est la mise en série des deux souplesses, \(1/(L^3/12 E I + L/G A_s)\). \(\Phi = 0\) redonne Euler-Bernoulli terme à terme.

Parce que la base dépend du matériau, aucun espace ne peut la tabuler : l’espace déclare MODEL_EMBEDDED, une interpolation que la physique porte elle-même.

La masse est intégrée des mêmes fonctions de forme : \(\int \rho A\,N_w^{\mathsf T} N_w + \int \rho I\,N_\theta^{\mathsf T} N_\theta\), le second terme étant l’inertie de rotation de la section ; l’axial et la torsion gardent la masse linéaire \((\rho L/6)\,[[2, 1], [1, 2]]\). À \(\Phi = 0\), elle redonne la table classique \(\rho A L/420\cdot[156, 22L, 54, -13L, \dots]\).

Aux points de Gauss, la courbure varie linéairement et le cisaillement est constant sur une travée sans charge, comme l’imposent \(M' = V\) et \(V' = 0\). Les reconstruire demande le matériau, puisque \(\Phi\) en dépend.

Variables

La sortie porte \(\Phi\) (phi) à côté des efforts : le noyau des forces internes reçoit l’état, pas le matériau, et son \(B\) dépend de \(\Phi\).

Ce qu’elle produit

Exemple

La console de la page Euler-Bernoulli, avec une aire de cisaillement : la flèche gagne \(F\,L/(G A_s)\).

poutre = pyrucast.model.timoshenko(
    pyrucast.FiniteElementSpace(ligne, interpolation="MODEL_EMBEDDED")
)
modele = poutre
for v in ("u_x", "u_y", "r_z"):
    modele = modele | pyrucast.model.dirichlet(encastrement, poutre, v)
materiaux = pyrucast.element_field.material_field(
    modele, [("E", 200e3), ("A", 1e-2), ("I", 1e-4), ("G", 80e3), ("A_s", 5e-3)]
)
u = pyrucast.matrix.stiffness(modele, materiaux).solve(force)
print(round(u.value(b, "u_y"), 12))  # -1.383333333333 = F·L³/(3EI) + F·L/(G·A_s)

Côté Rust

Timoshenko (src/models/timoshenko.rs) partage avec Bernoulli le \(B\) de src/models/beam.rs. Le test compare une console à \(P L^3/(3 E I) + P L/(G A_s)\) :

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::model;
use pyrucast::ops::solver::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))?);
        model::dirichlet(&imposed, target, var, 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::model;
use pyrucast::ops::solver::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))?);
        model::dirichlet(&imposed, target, var, 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(())
}

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::model;
use pyrucast::ops::solver::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))?);
        model::dirichlet(&imposed, target, var, 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(())
}