Barre / treillis

Élément SEG2 à 2 nœuds transmettant uniquement l’effort axial (treillis). Fonctionne à l’identique en 1-D, 2-D et 3-D.

Équations continues résolues

La barre suit la loi axiale 1-D et son équilibre le long de l’axe s :

\[ N = E\,A\,\varepsilon, \qquad \varepsilon = \frac{du}{ds}, \qquad \frac{dN}{ds} + f = 0, \]

où N est l’effort normal, ε la déformation axiale (dérivée du déplacement le long de la barre) et f la charge axiale répartie. L’orientation est déduite des coordonnées des nœuds via le cosinus directeur c = (x_B − x_A)/L.

Forme discrétisée

Avec l’interpolation linéaire SEG2, la déformation est constante par élément. Projetée sur les directions physiques par c, la rigidité élémentaire globale (en d dimensions) s’écrit

\[ K_e = \frac{E\,A}{L} \begin{bmatrix} c\,c^\top & -c\,c^\top \\ -c\,c^\top & c\,c^\top \end{bmatrix}, \]

écrite aux positions (NodeId_i, f_a) × (NodeId_j, u_b). En 1-D, c = 1 et l’on retrouve (EA/L)[[1,−1],[−1,1]].

Variables et matériau

• primal : u_x, u_y(, u_z) — dual : f_x, f_y(, f_z).
• matériau : E (module d’Young), A (section) ; rho facultatif (masse).
• comportement (COMP) : effort axial N = E·A·(cᵀ ε c), à partir de la déformation ε (op deformation).

⚠️ Une barre n’a aucune raideur transversale : pour un système bien posé il faut bloquer les DOFs transverses (treillis triangulé, appuis), sinon la matrice est singulière.

Mise en donnée (Rust, testé)

Barre horizontale encastrée à gauche, appuyée transversalement à droite, force axiale F au bout ⇒ u_x = F·L/(E·A). Le code est le test d’intégration tests/truss.rs, exécuté à chaque cargo test :

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::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 truss_bar_recovers_axial_elongation() -> Result<()> {
    const E: f64 = 210.0e9; // Young's modulus (Pa)
    const A: f64 = 1.0e-4; // section area (m²)
    const L: f64 = 2.0; // length (m)
    const F: f64 = 1000.0; // axial force at the right end (N)

    // ── Mesh: one horizontal SEG2 bar ──────────────────────────────────────
    let coords = Handle::new(Coords::new(2)?);
    let n0 = Node::create_in(coords.clone(), &[0.0, 0.0])?;
    let n1 = Node::create_in(coords.clone(), &[L, 0.0])?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    mesh.add_cell(&[n0.id(), n1.id()])?;
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : barre + appuis (Dirichlet homogènes) ──────────────────────
    // Homogeneous (u = 0) BCs: the imposed value defaults to 0, so we only need
    // to introduce the constraint. A bar has no transverse stiffness, hence
    // `u_y` is clamped at both nodes to make the system well-posed.
    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::truss(&fes)?;
    model = model.union(&clamp(&model, &n0, "u_x")?)?;
    model = model.union(&clamp(&model, &n0, "u_y")?)?;
    model = model.union(&clamp(&model, &n1, "u_y")?)?;

    // ── Matériau E, A (Dirichlet ignoré automatiquement) ───────────────────
    let materials = pyrucast::ops::element_field::material_field(&model, &[("E", E), ("A", A)])?;

    // ── Loading: axial force F at the right node ───────────────────────────
    let mut load_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    load_sm.add_cell(&[n1.id()])?;
    let load_sm = Handle::new(load_sm);
    let mut rhs = SubNodeField::from_poi1(&load_sm, vec!["f_x".into()])?;
    rhs.set_value(n1.id(), "f_x", F)?;
    let rhs = NodeField::from_sub(rhs);

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

    // ── Comparaison à l'analytique : u_x = F·L / (E·A) ─────────────────────
    let expected = F * L / (E * A);
    let ux = solution.value(n1.id(), "u_x")?;
    assert!(
        (ux - expected).abs() < 1e-10 * expected,
        "u_x = {ux}, attendu {expected}"
    );
    // The left node is clamped, the right end does not move transversally.
    assert!(solution.value(n0.id(), "u_x")?.abs() < 1e-18);
    assert!(solution.value(n1.id(), "u_y")?.abs() < 1e-18);

    Ok(())
}

Exemple Python

"""Bar / truss — a bar in tension, compared with the analytical solution.

Physique
--------
A 2-node `SEG2` element transmitting the axial force only. Law: `N = E·A·ε`
with `ε = du/ds` (axial strain along the bar). Global stiffness
`K_e = (E·A/L)·[[c⊗c, -c⊗c], [-c⊗c, c⊗c]]`, where `c` is the direction cosine
(derived from the nodes' coordinates) — works in 1-D/2-D/3-D.

Problème
--------
A horizontal bar of length `L`, clamped on the left (`u_x = u_y = 0`),
transversally supported on the right (`u_y = 0`), axial force `F` on the right.
A bar having no transverse stiffness, `u_y` is blocked at both nodes.
Solution analytique : `u_x = F·L / (E·A)`.

Lancement ::

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

import pyrucast

E, A, L, F = 210.0e9, 1.0e-4, 2.0, 1000.0


def _clamp(target, node, var):
    """Homogeneous Dirichlet (u = 0) on `var` at node `node`."""
    imposed = pyrucast.mesh.poi1_from_nodes([node])
    multiplier = pyrucast.mesh.barycenter(imposed)
    return pyrucast.model.dirichlet(target, var, imposed, multiplier)


def main() -> None:
    c = pyrucast.Coords(2)
    n0 = c.add_node([0.0, 0.0])
    n1 = c.add_node([L, 0.0])
    mesh = pyrucast.mesh.line(n0, n1, 1)  # un seul SEG2 (mailleur `line`)
    fes = pyrucast.FiniteElementSpace(mesh)

    model = pyrucast.model.truss(fes)
    model = model | _clamp(model, n0, "u_x")
    model = model | _clamp(model, n0, "u_y")
    model = model | _clamp(model, n1, "u_y")  # no transverse stiffness

    materials = pyrucast.element_field.material_field(model, [("E", E), ("A", A)])

    load = pyrucast.mesh.poi1_from_nodes([n1])
    rhs = pyrucast.NodeField(load, ["f_x"])
    rhs[0].set_value(n1, "f_x", F)

    solution = pyrucast.solver.solve(pyrucast.matrix.stiffness(model, materials), rhs)

    ux = solution.value(n1, "u_x")
    expected = F * L / (E * A)
    print(f"u_x (bout) = {ux:.6e}   (analytique F·L/E·A = {expected:.6e})")
    assert abs(ux - expected) < 1e-10 * expected
    print("OK : élongation axiale conforme à F·L/(E·A).")


if __name__ == "__main__":
    main()

Compléments

Masse & rigidité géométrique

En plus de la rigidité, la barre assemble :

• masse consistante M = (ρAL/6)[[2,1],[1,2]] sur chaque composante de translation (rho composante matériau optionnelle) — pyrucast.matrix.mass ;
• rigidité géométrique K_g = (N/L)·(I − c⊗c) transverse, sous l’effort axial N (sortie n du comportement) — pyrucast.matrix.geometric.