Poutre d’Euler-Bernoulli

Introduction

La théorie classique des poutres : les sections planes restent planes et normales à l’axe déformé, si bien que la rotation de section est la pente, θ = w', et qu’il n’y a aucun cisaillement transverse. C’est là toute la différence avec Timoshenko, le portique 2D et le cadre 3D, qui conservent une souplesse de cisaillement.

Trois configurations partagent la même théorie, et la dimension du maillage les départage — il n’y a rien à choisir :

Ce fut un argument, qui ne pouvait prendre que la valeur correspondant au maillage : toute autre était refusée. Un argument à valeur unique ne transporte aucune information — il n’offre qu’un moyen de se contredire — donc la configuration se déduit. Les noms de DDL obtenus se relisent par model.primal_vars().

Équations continues résolues

La cinématique est celle d’une section rigide qui tourne avec la pente :

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

D’où une déformation axiale affine dans la section, et aucune distorsion :

\[ \varepsilon_{xx} = u’ - y\,w’’ = \varepsilon_0 - y\,\chi, \qquad \chi = w’’, \qquad \gamma_{xy} \equiv 0 . \]

\( \chi \) est la courbure. En intégrant \( \sigma_{xx} = E\varepsilon_{xx} \) sur la section, les efforts généralisés se découplent — c’est le choix de l’axe neutre, \( \int_A y\,dA = 0 \) — et la loi de section s’écrit

\[ N = EA\,\varepsilon_0, \qquad M = EI\,\chi, \qquad I = \int_A y^2\,dA . \]

L’équilibre local intégré sur la section donne alors les deux équations classiques, découplées elles aussi :

\[ (EA\,u’)’ + n = 0, \qquad (EI\,w’‘)’’ = q . \]

La seconde est du quatrième ordre — et c’est toute la différence avec Timoshenko, qui en fait deux du second ordre en gardant \( \theta \) indépendant de \( w’ \). En 3-D s’y ajoutent la flexion selon le second axe principal (\( M_z = EI_z\,w_y’’ \)) et la torsion de Saint-Venant \( M_t = GJ\,\varphi’ \), qui ne se couple à rien pour une section symétrique.

Pourquoi une physique à part, et non une aire de cisaillement infinie

On pourrait atteindre Bernoulli en faisant tendre G·A_s → ∞ dans un élément de Timoshenko, et le résultat serait juste en arithmétique exacte. En virgule flottante il ne l’est pas : le terme de cisaillement domine alors la raideur de plusieurs ordres de grandeur et la réponse en flexion s’y noie — le blocage en cisaillement classique, atteint par l’autre bout.

Écrire la théorie directement supprime la question, et supprime au passage deux constantes matériau (G, A_s) qu’un modèle de Bernoulli n’a aucune raison de demander. Réclamer une constante qu’une théorie n’utilise pas, c’est inviter la mauvaise.

L’élément

L’équation étant du quatrième ordre, sa forme faible demande une interpolation \( C^1 \) : le déplacement et sa pente doivent être continus d’un élément au suivant. C’est exactement ce que fournit la famille Hermite3, qui prend pour degrés de liberté la flèche et la pente à chaque extrémité — deux fonctions de forme par nœud au lieu d’une.

Ce n’est donc pas une interpolation de Lagrange, et le modèle exige un espace HERMITE3 :

fes = pyrucast.FiniteElementSpace(maillage, interpolation="HERMITE3")
poutre = pyrucast.model.bernoulli(fes)  # 1-D ⇒ flexion pure

Un espace de Lagrange porterait une flèche linéaire, de courbure identiquement nulle : le modèle le refuse plutôt que d’assembler une raideur qui ne correspondrait pas à la base déclarée.

La courbure est alors linéaire sur l’élément, donc l’espace d’approximation contient exactement la solution d’une travée chargée seulement à ses extrémités : l’élément est exact aux nœuds pour toute charge laissant la travée libre d’efforts répartis — c’est pourquoi un élément par barre suffit pour un portique. La raideur \( K_b = \int_0^L EI\,N’‘^\top N’’\,dx \) est intégrée depuis cette base, et vaut exactement la forme fermée classique :

\[ K_b = \frac{EI}{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 \text{DDL } [\,w_A,\ \theta_A,\ w_B,\ \theta_B\,]. \]

L’intégration plutôt que la forme fermée est un choix : elle laisse une seule source de vérité, et rend l’interpolation déclarée porteuse. Une base fausse produirait désormais une raideur fausse, que les tests de poutre attraperaient ; avec une matrice écrite en dur, elle aurait pu être n’importe quoi. La forme fermée reste, comme oracle de test — c’est elle que tests/hermite.rs compare à l’intégrale, à la précision machine.

L’effort axial y est ajouté par le terme de barre \( EA/L \), qui ne s’y couple pas.

Le repère local 3-D est déduit automatiquement d’une référence globale Z (globale Y pour une barre quasi verticale), comme pour le cadre 3D : aucune donnée d’orientation à fournir, ce qui convient aux sections symétriques.

Variables et matériau

Le comportement (COMP) rend les efforts de section — M en 1-D, N, M en plan, N, M_y, M_z, T dans l’espace — par une loi linéaire, comme tout élément structural.

Elle n’en rend aucun de cisaillement, et la liste s’arrête donc plus tôt que celle de Timoshenko. Ce n’est pas une omission : Euler-Bernoulli est Φ = 0, donc les lignes de cisaillement de son B sont nulles, et il n’y a rien à leur apparier. Les deux absences sont un seul énoncé — c’est aussi pourquoi cette théorie n’a pas à porter Φ dans son état, là où l’autre le doit.

Les forces internes intègrent le transposé du même B (models::beam::b_into) que la rigidité : ∫ Bᵀσ vaut K·u exactement, la loi étant linéaire. C’est ce que mesure tests/internal_forces.rs, dans les trois configurations.

Mise en donnée (Rust, testé)

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::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::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))?);
        let multiplier = mesh::barycenter(&imposed)?;
        model = model.union(&model::dirichlet(
            &model,
            var,
            &imposed,
            &multiplier,
            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(())
}

Exemple Python

model = pyrucast.model.bernoulli(fes)  # Coords 2-D ⇒ portique plan
materials = pyrucast.element_field.material_field(
    model, [("E", 210_000.0), ("A", 1e-2), ("I", 1e-4)]
)
k = pyrucast.matrix.stiffness(model, materials)

Compléments

Masse et rigidité géométrique. La masse cohérente vient du même bloc que celui de Timoshenko, pris à Φ = 0 : Bernoulli n’apporte donc aucune dérivation propre, il est le bout sans cisaillement d’une seule. Ce que ce bloc rend à Φ = 0 est la table classique ρAL/420·[156, 22L, 54, −13L ; …], ce qu’un test affirme.

La rigidité géométrique demande un effort axial pour raidir la barre : la configuration 1-D, en flexion pure, n’en déclare donc aucune — elle ne contribue rien plutôt que d’erroner, un modèle qui ne déclare pas un genre étant simplement ignoré de l’assembleur.

Ce que valent les tests. Un élément de poutre gagne sa place en étant exact aux nœuds : les tests comparent donc aux formules du cours à la précision machine (1e-12), et non à une tolérance de discrétisation. Console sous charge en bout (PL³/3EI), sous moment en bout (ML²/2EI — un cas qu’un signe faux dans la matrice d’Hermite raterait tout en passant le premier), traction axiale découplée de la flexion, et torsion TL/GJ.

Quand préférer Timoshenko. Un dernier test mesure ce qui sépare les deux théories : une poutre élancée donne le même résultat aux deux (à 1 % près), une poutre trapue fléchit nettement plus avec le cisaillement. Bernoulli est exactement la théorie qui dit qu’elle ne le fait pas — à utiliser tant que l’élancement le permet, et à quitter sinon.

Ce test-là maille les deux poutres. Bernoulli est exact avec un seul élément, mais l’interpolation linéaire de Timoshenko ne l’est pas : comparer les deux théories sur un élément unique mesurerait le maillage, pas la physique.