Poutre d’Euler-Bernoulli

En bref

La théorie classique des poutres : les sections restent planes et normales à l’axe déformé, si bien que la rotation de la section est la pente, et qu’il n’y a aucun cisaillement transverse.

La configuration se lit sur la dimension des Coords :

Théorie

La section rigide tourne avec la pente :

\[ u_x(x, y) = u(x) - y\,w'(x), \qquad u_y = w(x), \qquad \theta = w' , \]

d’où une déformation axiale affine dans la section, \(\varepsilon_{xx} = \varepsilon_0 - y\,\chi\), de courbure \(\chi = w''\), et aucune distorsion. Intégrée sur la section, la loi de Hooke donne les efforts de section,

\[ N = E A\,\varepsilon_0, \qquad M = E I\,\chi , \]

et l’équilibre deux équations découplées :

\[ (E A\,u')' + n = 0, \qquad (E I\,w'')'' = q . \]

La seconde est du quatrième ordre : sa forme faible demande une interpolation \(C^1\), le déplacement et sa pente continus d’un élément au suivant. C’est ce que donne l’interpolation HERMITE3, qui prend la flèche et la pente à chaque extrémité. Un espace de Lagrange est refusé : sa flèche linéaire aurait une courbure nulle.

La courbure est linéaire sur l’élément : il est exact aux nœuds pour toute travée chargée à ses seules extrémités, et un élément par barre suffit. La raideur de flexion \(\int_0^L E I\,N''^{\mathsf T} N''\) est intégrée depuis la base et vaut la forme fermée classique,

\[ K_b = \frac{E I}{L^3} \begin{bmatrix} 12 & 6L & -12 & 6L \\ 6L & 4L^2 & -6L & 2L^2 \\ -12 & -6L & 12 & -6L \\ 6L & 2L^2 & -6L & 4L^2 \end{bmatrix}, \qquad [\,w_A,\ \theta_A,\ w_B,\ \theta_B\,]. \]

L’effort axial s’y ajoute par le terme de barre \(E A/L\), sans couplage. En 3-D, la flexion selon le second axe principal (\(M_z = E I_z\,\chi_z\)) et la torsion de Saint-Venant (\(T = G J\,\varphi'\)) s’y ajoutent ; le repère local est déduit d’une référence globale Z (ou Y pour une barre presque verticale), sans donnée d’orientation, ce qui convient aux sections symétriques.

Bernoulli est la limite sans cisaillement de Timoshenko. Une physique à part évite d’y arriver par \(G A_s \to \infty\), où le terme de cisaillement noierait la flexion en virgule flottante, et évite de demander G et A_s, que cette théorie n’utilise pas.

Variables

L’entrée vient de pyrucast.element_field.beam_deformation, ou de strain.

Ce qu’elle produit

Exemple

Une console de longueur 2 dans le plan, encastrée, chargée en bout. La flèche est exacte, \(F\,L^3/(3 E I)\) :

# Une console de longueur 2 dans le plan, encastrée à gauche, F = −10 en y au bout.
c = pyrucast.Coords(2)
a = c.add_node([0.0, 0.0])
b = c.add_node([2.0, 0.0])
ligne = pyrucast.mesh.line(a, b, 2)
poutre = pyrucast.model.bernoulli(
    pyrucast.FiniteElementSpace(ligne, interpolation="HERMITE3")
)
print(poutre[0].primal_vars())  # ['u_x', 'u_y', 'r_z'] : un portique plan
encastrement = pyrucast.mesh.poi1_from_nodes([a])
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)]
)

force = pyrucast.NodeField(ligne, ["f_y"])
force.unit().set_value(b, "f_y", -10.0)
u = pyrucast.matrix.stiffness(modele, materiaux).solve(force)
print(round(u.value(b, "u_y"), 12))  # -1.333333333333 = F·L³/(3EI), exacte

Côté Rust

Bernoulli (src/models/bernoulli.rs) refuse à la construction un espace qui n’est pas HERMITE3. Son \(B\) est celui de src/models/beam.rs, partagé avec Timoshenko pour les forces internes. tests/hermite.rs compare la raideur intégrée à la forme fermée, à la précision machine ; le test ci-dessous compare la console sous charge et sous moment en bout, la traction et la torsion aux formules du cours, à 1e-12 :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Interpolation, Node};
use pyrucast::containers::field::SubField;
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;

const E: f64 = 210_000.0;
const I: f64 = 1.0e-4;
const L: f64 = 2.0;

#[test]
fn a_cantilever_under_a_tip_load_matches_its_closed_form() -> Result<()> {
    const P: f64 = 50.0;
    let coords = Handle::new(Coords::new(1)?);
    let a = Node::create_in(coords.clone(), &[0.0])?;
    let b = Node::create_in(coords.clone(), &[L])?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    mesh.add_cell(&[a.id(), b.id()])?;
    let fes = FiniteElementSpace::new(&mesh, Interpolation::Hermite3)?;

    // Clamped at A: both the deflection and the rotation are held.
    let mut model = model::bernoulli(&fes)?;
    for var in ["w", "theta"] {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(&a))?);
        model = model.union(&model::dirichlet(
            &imposed,
            &model,
            var,
            Default::default(),
        )?)?;
    }
    let materials = pyrucast::ops::element_field::material_field(&model, &[("E", E), ("I", I)])?;

    // A point load at the free end.
    let load_sm = Handle::new(SubMesh::poi1_from_nodes(std::slice::from_ref(&b))?);
    let mut rhs = SubNodeField::from_poi1(&load_sm, vec!["f_w".into()])?;
    rhs.set_value(b.id(), "f_w", P)?;

    let k = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&k, &NodeField::from_sub(rhs))?;

    // Nodally exact: w = PL³/3EI, θ = PL²/2EI, to machine precision.
    let w = solution.value(b.id(), "w")?;
    let theta = solution.value(b.id(), "theta")?;
    let w_exact = P * L.powi(3) / (3.0 * E * I);
    let theta_exact = P * L * L / (2.0 * E * I);
    assert!(
        (w - w_exact).abs() < 1e-12 * w_exact,
        "w = {w}, exact {w_exact}"
    );
    assert!(
        (theta - theta_exact).abs() < 1e-12 * theta_exact,
        "θ = {theta}, exact {theta_exact}"
    );
    Ok(())
}