Élasticité linéaire
Continuum en petites déformations : 2-D (TRI3 / QUA4) ou 3-D
(TET4 / HEX8), les cas 2-D couvrant aussi l’axisymétrie (solide de
révolution maillé dans son plan méridien).
Équations continues résolues
Sur le domaine \( \Omega \), avec \( b \) les efforts volumiques :
\[ \underbrace{\nabla\cdot\sigma + b = 0}{\text{équilibre}}, \qquad \underbrace{\sigma = \mathbb{D} : \varepsilon}{\text{loi de Hooke}}, \qquad \underbrace{\varepsilon = \tfrac12(\nabla u + \nabla u^\top)}_{\text{cinématique}}. \]
La forme faible (multiplication par un déplacement virtuel \( v \), intégration par parties) s’écrit : trouver \( u \) tel que pour tout \( v \),
\[ \int_\Omega \varepsilon(v) : \mathbb{D} : \varepsilon(u)\,d\Omega = \int_\Omega v\cdot b\,d\Omega + \int_{\Gamma_N} v\cdot t\,d\Gamma, \]
où \( t \) est la traction imposée sur le bord de Neumann \( \Gamma_N \).
Forme discrétisée
En convention de Voigt (déformation ingénieur \( \gamma = 2\varepsilon \)),
le champ discret \( u_h = \sum_i N_i u_i \) donne \( \varepsilon = B\,u_e \),
avec la matrice déformation-déplacement \( B \) bâtie des dérivées
physiques \( \partial N_i/\partial x_a \) (voir
dn_dx). En 2-D
(\( \varepsilon = [\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy}]^\top \)),
le bloc du nœud \( i \) est
\[ B_i = \begin{bmatrix} \partial_x N_i & 0 \\ 0 & \partial_y N_i \\ \partial_y N_i & \partial_x N_i \end{bmatrix}, \]
et en 3-D (\( \varepsilon = [\varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz}, \gamma_{yz}, \gamma_{xz}, \gamma_{xy}]^\top \)),
\[ B_i = \begin{bmatrix} \partial_x N_i & 0 & 0 \\ 0 & \partial_y N_i & 0 \\ 0 & 0 & \partial_z N_i \\ 0 & \partial_z N_i & \partial_y N_i \\ \partial_z N_i & 0 & \partial_x N_i \\ \partial_y N_i & \partial_x N_i & 0 \end{bmatrix}. \]
La rigidité élémentaire est alors, intégrée par quadrature de Gauss,
\[ K_e = \int_{\Omega_e} B^\top D\, B\, d\Omega \;\approx\; \sum_g B(\xi_g)^\top D\, B(\xi_g)\,|J(\xi_g)|\,w_g, \]
écrite aux positions (NodeId_i, f_a) × (NodeId_j, u_b) (ordre des DOFs
nœud-majeur). Le second membre nodal cohérent d’une traction de bord est
\( f_i = \int_{\Gamma_N} N_i\,t\,d\Gamma \) (opérateur
flux).
Matrice constitutive D
Le modèle fixe \( D \) (isotrope, module d’Young \( E \), coefficient de Poisson \( \nu \)) :
plane_stress(contraintes planes, \( \sigma_{zz}=0 \)), avec \( c = \dfrac{E}{1-\nu^2} \) :
\[ D = c\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \tfrac{1-\nu}{2} \end{bmatrix}; \]
plane_strain(déformations planes, \( \varepsilon_{zz}=0 \), \( \sigma_{zz}\neq 0 \)), avec \( c = \dfrac{E}{(1+\nu)(1-2\nu)} \) :
\[ D = c\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \tfrac{1-2\nu}{2} \end{bmatrix}; \]
full_3d(3-D), même \( c \), avec le module de cisaillement \( G = c\,\tfrac{1-2\nu}{2} \) :
\[ D = \begin{bmatrix} c(1-\nu) & c\nu & c\nu & & & \\ c\nu & c(1-\nu) & c\nu & & & \\ c\nu & c\nu & c(1-\nu) & & & \\ & & & G & & \\ & & & & G & \\ & & & & & G \end{bmatrix} \quad (\text{ordre } [xx, yy, zz, yz, xz, xy]). \]
axisymmetric(solide de révolution), même \( c \) : les trois directions normales \( r, z, \theta \) étant orthogonales, le bloc normal est l’isotrope 3×3 et \( rz \) est le seul cisaillement,
\[ D = c\begin{bmatrix} 1-\nu & \nu & \nu & \\ \nu & 1-\nu & \nu & \\ \nu & \nu & 1-\nu & \\ & & & \tfrac{1-2\nu}{2} \end{bmatrix} \quad (\text{ordre } [rr, zz, \theta\theta, rz]). \]
Axisymétrie
Un solide de révolution se maille dans son plan méridien \( (r, z) \) sur
des Coords déclarées axisymétriques (\( x = r \ge 0 \), \( y = z \)) — voir
Coordonnées. Deux choses changent, et elles
ont deux origines distinctes :
- la mesure d’intégration, portée par la géométrie :
\( d\Omega = 2\pi r\,|J|\,d\xi \). Elle vaut pour toutes les intégrales
— rigidité, masse, conductivité, flux réparti, volumes, forces internes, y
compris sur les sous-maillages de bord
SEG2, dont \( \int 2\pi r\,N \) donne directement l’effort sur l’anneau. Rien à écrire : c’estCellGeom::det_j_wqui l’applique, en un seul point ; - la déformation orthoradiale, portée par le modèle : \( \varepsilon_{\theta\theta} = u_r / r \), que le gradient méridien ne peut pas exprimer. Elle ajoute une quatrième composante de Voigt et une ligne à \( B \) :
\[ B_i = \begin{bmatrix} \partial_r N_i & 0 \\ 0 & \partial_z N_i \\ N_i / r & 0 \\ \partial_z N_i & \partial_r N_i \end{bmatrix} \quad (\varepsilon = [\varepsilon_{rr}, \varepsilon_{zz}, \varepsilon_{\theta\theta}, \gamma_{rz}]^\top). \]
Les points de Gauss étant intérieurs à la maille, \( r > 0 \) même pour un élément qui touche l’axe : le terme \( N_i/r \) reste fini, sans traitement particulier de l’axe.
Nommage (convention Cast3M) : les composantes s’appellent sigma_xx,
sigma_yy, sigma_zz, sigma_xy et eps_xx, eps_yy, eps_zz, eps_xy,
où zz désigne l’orthoradial \( \theta\theta \) — le plan méridien
n’occupant que xx, yy et xy, il n’y a pas de collision.
Le modèle et le repère doivent s’accorder dans les deux sens : une géométrie
de révolution refuse plane_stress / plane_strain, et axisymmetric refuse
une géométrie cartésienne. Sans cela on mélangerait silencieusement une loi plane
avec la mesure \( 2\pi r \).
La thermique n’a rien de spécifique à faire : le flux \( q = -k\nabla T \) est déjà purement méridien, et le facteur \( 2\pi r \) suffit à produire le profil logarithmique d’un cylindre creux. La plasticité et Mazars supportent l’axisymétrie, leur état interne étant déjà stocké en 3-D complet. En revanche barre et portique la refusent : un segment du plan méridien engendre une coque de révolution, que leurs noyaux ne modélisent pas.
Un maillage de bord (SEG2 en 2-D) est par ailleurs refusé comme domaine
par les trois physiques de milieu continu : B y serait bâti sur le gradient
tangent et \( B^\top D B \) serait déficient en rang dans la direction
normale. Un bord porte des charges (flux, convection), il n’est pas un massif.
Validation : tests/axisymmetric.rs (Lamé, patch test de dilatation uniforme,
\( \int B^\top\sigma = K u \), volume et masse de révolution, conduction
logarithmique) et tests/python/test_axisymmetric.py.
Convergence sur Lamé
La solution de Lamé \( u_r = c_1 r + c_2/r \) comporte un terme rationnel : aucune base de Lagrange, de quelque degré que ce soit, ne la reproduit exactement. Les éléments quadratiques gagnent un ordre, pas l’exactitude. La solution ne dépendant que de \( r \), le problème discret est une EDO 1-D et les valeurs nodales sont superconvergentes en \( O(h^{2p}) \) :
nr | Q1 (QUA4) | ordre | Q2 (QUA8) | ordre |
|---|---|---|---|---|
| 5 | 6,5e-3 | — | 2,6e-5 | — |
| 10 | 1,6e-3 | 1,97 | 1,7e-6 | 3,94 |
| 20 | 4,1e-4 | 1,99 | 1,1e-7 | 3,99 |
| 40 | 1,0e-4 | 2,00 | 6,8e-9 | 4,00 |
(erreur relative maximale sur \( u_r \)). Les contraintes, une dérivée plus bas, passent de \( O(h) \) à \( O(h^2) \).
Le cas exact existe néanmoins : lorsque \( c_2 = 0 \) — dilatation uniforme \( u_r = c\,r \) — l’état de déformation est constant et même Q1 le reproduit à la précision machine (c’est le patch test de la suite de validation).
Matrice de masse
Pour la dynamique, la masse consistante (composante matériau rho) est
\[ M_e = \int_{\Omega_e} \rho\,N^\top N\, d\Omega \;\approx\; \sum_g \rho\,N(\xi_g)^\top N(\xi_g)\,|J(\xi_g)|\,w_g, \]
où \( N \) place \( N_i \) sur chaque composante de translation — assemblée
par assemble.mass, et concentrable en diagonale
par lump.
Variables et matériau
- primal :
u_x, u_y(, u_z)— dual :f_x, f_y(, f_z). - matériau :
E(Young),nu(Poisson) ; facultatifalpha(dilatation thermique, cf. thermomécanique),rho(masse) — accepté par le champ matériau mais jamais exigé pour un assemblage purement élastique. - comportement (
COMP) :σ = D ε(convention tenseur → ingénieurγ = 2ε), à partir de la déformationε(opdeformation). - modèles :
plane_stress,plane_strain,axisymmetric(2-D) etfull_3d(3-D).
Mise en donnée (Rust, testé)
Carré unité en contraintes planes : appuis u_x = 0 (gauche) et u_y = 0
(bas), traction S sur le bord droit appliquée en charges nodales cohérentes
par l’opérateur flux (composante f_x).
Solution exacte u_x = (S/E)·x, u_y = −(ν S/E)·y. Code = test
tests/elasticity.rs (le fichier contient aussi un test 3-D sur un cube
HEX8) :
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::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;
#[test]
fn elasticity_unit_square_uniaxial_tension() -> Result<()> {
const E: f64 = 210.0; // Young's modulus
const NU: f64 = 0.3; // Poisson's ratio
const S: f64 = 2.0; // traction on the right edge
const N: usize = 2; // N×N QUA4 grid
let h = 1.0 / N as f64;
// ── QUA4 mesh on [0,1]² ────────────────────────────────────────────────
let coords = Handle::new(Coords::new(2)?);
let idx = |i: usize, j: usize| j * (N + 1) + i;
let mut grid: Vec<Node> = Vec::new();
for j in 0..=N {
for i in 0..=N {
grid.push(Node::create_in(
coords.clone(),
&[i as f64 * h, j as f64 * h],
)?);
}
}
let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
for j in 0..N {
for i in 0..N {
mesh.add_cell(&[
grid[idx(i, j)].id(),
grid[idx(i + 1, j)].id(),
grid[idx(i + 1, j + 1)].id(),
grid[idx(i, j + 1)].id(),
])?;
}
}
let fes = FiniteElementSpace::lagrange1(&mesh)?;
// ── Modèle : élasticité plane stress + appuis (rollers) ────────────────
// The support receives the model it constrains: the dual is read there, and
// contraints s'y vérifient.
let roller = |target: &Model, nodes: &[Node], var: &str| -> Result<Model> {
let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(nodes)?);
let multiplier = mesh::barycenter(&imposed)?;
model::dirichlet(target, var, &imposed, &multiplier, Default::default())
};
let left: Vec<Node> = (0..=N).map(|j| grid[idx(0, j)].clone()).collect();
let bottom: Vec<Node> = (0..=N).map(|i| grid[idx(i, 0)].clone()).collect();
let mut model = model::elasticity(&fes, Kinematics::PlaneStress)?;
model = model.union(&roller(&model, &left, "u_x")?)?;
model = model.union(&roller(&model, &bottom, "u_y")?)?;
// ── Loading: traction S on the right edge (consistent nodal loads, on the
// f_x component). The load is a sub-model: it joins the model, its density
// the material. ──────────────────────────────────────────────────────
let mut right_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
for j in 0..N {
right_edge.add_cell(&[grid[idx(N, j)].id(), grid[idx(N, j + 1)].id()])?;
}
let right_fes = FiniteElementSpace::lagrange1(&right_edge)?;
model = model.union(&model::flux(&right_fes, &model, "f_x".into())?)?;
let materials = pyrucast::ops::element_field::material_field(
&model,
&[("E", E), ("nu", NU), ("phi_f_x", S)],
)?;
let rhs = pyrucast::ops::node_field::external_forces(&model, &materials)?;
// ── Assemblage + résolution ────────────────────────────────────────────
let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
let solution = solve(&stiffness, &rhs)?;
// ── Compared with the analytical u_x = (S/E)·x, u_y = −(ν S/E)·y ───────
let tol = 1e-10;
for j in 0..=N {
for i in 0..=N {
let (x, y) = (i as f64 * h, j as f64 * h);
let ux = solution.value(grid[idx(i, j)].id(), "u_x")?;
let uy = solution.value(grid[idx(i, j)].id(), "u_y")?;
assert!((ux - S / E * x).abs() < tol, "u_x({x},{y}) = {ux}");
assert!((uy + NU * S / E * y).abs() < tol, "u_y({x},{y}) = {uy}");
}
}
Ok(())
}Exemple Python
"""Élasticité linéaire — traction d'un carré (contraintes planes).
Physique
--------
Continuum en petites déformations : équilibre `∇·σ = 0`, loi de Hooke
`σ = D : ε`, kinematics `ε = ½(∇u + ∇uᵀ)`. The stiffness is
`K = ∫ Bᵀ D B dΩ` (B : matrice déformation-déplacement en Voigt, D : matrice
constitutive isotrope, ici en contraintes planes).
Problème
--------
Carré unité, appuis `u_x = 0` (bord gauche) et `u_y = 0` (bord bas), traction
`S` on the right edge applied as consistent nodal loads by the `flux` operator
(on the `f_x` component). Exact (uniaxial) solution:
`u_x = (S/E)·x`, `u_y = -(ν·S/E)·y`.
Lancement ::
maturin develop --features extension-module
python examples/elasticity.py
"""
import pyrucast
E, NU, S, N = 210.0, 0.3, 2.0, 2
def _clamp(target, nodes, var):
imposed = pyrucast.mesh.poi1_from_nodes(nodes)
multiplier = pyrucast.mesh.barycenter(imposed)
return pyrucast.model.dirichlet(target, var, imposed, multiplier)
def main() -> None:
h = 1.0 / N
c = pyrucast.Coords(2)
def idx(i, j):
return j * (N + 1) + i
# An N×N grid of QUA4 by sweeping two SEG2 lines (`sweep`).
bottom = pyrucast.mesh.line(c.add_node([0.0, 0.0]), c.add_node([1.0, 0.0]), N)
top = pyrucast.mesh.line(c.add_node([0.0, 1.0]), c.add_node([1.0, 1.0]), N)
mesh = pyrucast.mesh.sweep(bottom, top, N)
# Nodes laid out by idx(i, j) (i along x, j along y) by reading the
# connectivité QUA4 : maille (cy, cx) = cy*N + cx, nœuds locaux 0..3.
grid = [None] * ((N + 1) * (N + 1))
for cy in range(N):
for cx in range(N):
cell = cy * N + cx
grid[idx(cx, cy)] = mesh.node(0, cell, 0)
grid[idx(cx + 1, cy)] = mesh.node(0, cell, 1)
grid[idx(cx + 1, cy + 1)] = mesh.node(0, cell, 2)
grid[idx(cx, cy + 1)] = mesh.node(0, cell, 3)
fes = pyrucast.FiniteElementSpace(mesh)
left = [grid[idx(0, j)] for j in range(N + 1)]
bottom = [grid[idx(i, 0)] for i in range(N + 1)]
model = pyrucast.model.elasticity(fes, "plane_stress")
model = model | _clamp(model, left, "u_x")
model = model | _clamp(model, bottom, "u_y")
# Traction S on the right edge → consistent nodal loads (the flux op).
right = pyrucast.Mesh(c, "SEG2")
for j in range(N):
right.unit().add_cell([grid[idx(N, j)], grid[idx(N, j + 1)]])
right_fes = pyrucast.FiniteElementSpace(right)
model = model | pyrucast.model.flux(right_fes, model, "f_x")
materials = pyrucast.element_field.material_field(
model, [("E", E), ("nu", NU), ("phi_f_x", S)]
)
rhs = pyrucast.node_field.external_forces(model, materials)
solution = pyrucast.solver.solve(pyrucast.matrix.stiffness(model, materials), rhs)
print(f"{'x':>5} {'y':>5} {'u_x':>12} {'u_y':>12}")
tol = 1e-10
for j in range(N + 1):
for i in range(N + 1):
x, y = i * h, j * h
ux = solution.value(grid[idx(i, j)], "u_x")
uy = solution.value(grid[idx(i, j)], "u_y")
print(f"{x:5.2f} {y:5.2f} {ux:12.6e} {uy:12.6e}")
assert abs(ux - S / E * x) < tol
assert abs(uy + NU * S / E * y) < tol
print("\nOK: uniaxial field matching u_x=(S/E)x, u_y=-(νS/E)y.")
if __name__ == "__main__":
main()
Compléments
Thermomécanique non couplée
Première brique de thermomécanique : une température imposée ΔT engendre une
déformation thermique de libre dilatation ε_th = α·(T − T_ref), d’où des
contraintes mécaniques — sans rétroaction de la mécanique sur le thermique.
En petites déformations, la rigidité K reste l’élastique ; le terme thermique
n’agit que sur le second membre et sur la contrainte réelle :
\[ \sigma = D : (\varepsilon(u) - \varepsilon_{th}), \qquad f_{th} = \int_\Omega B^\top D\, \varepsilon_{th}\, d\Omega. \]
Aucune physique nouvelle : on compose les briques existantes. alpha est fourni
au champ matériau (composante facultative) ; la température, portée aux points de
Gauss par interp_to_gauss, alimente
thermal_strain (EPTH) ; la charge thermique sort
de integrate_behavior + internal_forces (BSIG) ; enfin la contrainte réelle
se relit sur deformation(u) − ε_th.
Deux régimes sur une barre chauffée valident les fermetures analytiques : bord en
x encastré aux deux bouts ⇒ σ_xx = −E·α·ΔT ; appuis simples ⇒ dilatation libre
u = α·ΔT·(x, y) sans contrainte. Code = test tests/thermoelastic_bar.rs :
use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::field::Field;
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::models::tensor::Kinematics;
use pyrucast::ops::element_field::{deformation, interp_to_gauss, thermal_strain};
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::node_field::internal_forces;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;
#[test]
fn thermoelastic_constrained_bar_stress() -> Result<()> {
const E: f64 = 210_000.0;
const NU: f64 = 0.3;
const ALPHA: f64 = 1e-5;
const T_REF: f64 = 20.0;
const DT: f64 = 100.0;
const NX: usize = 4;
const NY: usize = 2;
const L: f64 = 4.0;
const H: f64 = 1.0;
let (hx, hy) = (L / NX as f64, H / NY as f64);
// ── QUA4 mesh on [0,L]×[0,H] ───────────────────────────────────────────
let coords = Handle::new(Coords::new(2)?);
let idx = |i: usize, j: usize| j * (NX + 1) + i;
let mut grid: Vec<Node> = Vec::new();
for j in 0..=NY {
for i in 0..=NX {
grid.push(Node::create_in(
coords.clone(),
&[i as f64 * hx, j as f64 * hy],
)?);
}
}
let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
for j in 0..NY {
for i in 0..NX {
mesh.add_cell(&[
grid[idx(i, j)].id(),
grid[idx(i + 1, j)].id(),
grid[idx(i + 1, j + 1)].id(),
grid[idx(i, j + 1)].id(),
])?;
}
}
let fes = FiniteElementSpace::lagrange1(&mesh)?;
// ── Modèle : élasticité + deux bords en x encastrés + appui u_y en bas ──
let clamp = |target: &Model, nodes: &[Node], var: &str| -> Result<Model> {
let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(nodes)?);
let multiplier = mesh::barycenter(&imposed)?;
model::dirichlet(target, var, &imposed, &multiplier, Default::default())
};
let left: Vec<Node> = (0..=NY).map(|j| grid[idx(0, j)].clone()).collect();
let right: Vec<Node> = (0..=NY).map(|j| grid[idx(NX, j)].clone()).collect();
let bottom: Vec<Node> = (0..=NX).map(|i| grid[idx(i, 0)].clone()).collect();
let mut model = model::elasticity(&fes, Kinematics::PlaneStress)?;
model = model.union(&clamp(&model, &left, "u_x")?)?;
model = model.union(&clamp(&model, &right, "u_x")?)?;
model = model.union(&clamp(&model, &bottom, "u_y")?)?;
// `alpha` supplied through the material field — an optional elastic component.
let materials = pyrucast::ops::element_field::material_field(
&model,
&[("E", E), ("nu", NU), ("alpha", ALPHA)],
)?;
// ── Imposed temperature T = T_ref + ΔT everywhere, carried at the Gauss points
let support = Handle::new(SubMesh::poi1_from_nodes(&grid)?);
let mut t_nodal = SubNodeField::from_poi1(&support, vec!["T".into()])?;
for n in &grid {
t_nodal.set_value(n.id(), "T", T_REF + DT)?;
}
let t_elem = interp_to_gauss(&NodeField::from_sub(t_nodal), &fes)?;
// ── Charge thermique f_th = ∫ Bᵀ D ε_th (BSIG de σ_th = D:ε_th) ─────────
let eps_th = thermal_strain(&t_elem, &materials, &fes, T_REF)?;
let sig_th =
pyrucast::ops::element_field::behavior::integrate(&model, &eps_th, None, &materials, None)?;
let f_th = internal_forces(&model, &sig_th, &NodeField::empty(), &materials)?;
// ── Assemblage + résolution ────────────────────────────────────────────
let solution = solve(
&pyrucast::ops::matrix::stiffness(&model, &materials)?,
&f_th,
)?;
// ── Déplacement propre (u_x, u_y) puis σ = D:(ε(u) − ε_th) ─────────────
let disp_support = Handle::new(SubMesh::poi1_from_nodes(&grid)?);
let mut disp = SubNodeField::from_poi1(&disp_support, vec!["u_x".into(), "u_y".into()])?;
for n in &grid {
disp.set_value(n.id(), "u_x", solution.value(n.id(), "u_x")?)?;
disp.set_value(n.id(), "u_y", solution.value(n.id(), "u_y")?)?;
}
let eps = deformation(&NodeField::from_sub(disp), &fes)?;
let eps_mech = eps.merge_field(&eps_th, |a, b| a - b)?;
let sigma = pyrucast::ops::element_field::behavior::integrate(
&model, &eps_mech, None, &materials, None,
)?;
// ── Vérification : σ_xx = −E·α·ΔT, σ_yy = 0 ────────────────────────────
let expected = -E * ALPHA * DT;
let tol = 1e-6 * expected.abs();
let sub = sigma.get(0)?.read();
for cell in 0..sub.cell_count() {
for g in 0..sub.gauss_count() {
assert!((sub.value(cell, g, "sigma_xx")? - expected).abs() < tol);
assert!(sub.value(cell, g, "sigma_yy")?.abs() < tol);
}
}
Ok(())
}"""Thermomécanique non couplée — barre chauffée (contraintes planes).
Physique
--------
An imposed temperature ΔT generates a free thermal strain
`ε_th = α·(T − T_ref)`. En petites déformations la rigidité reste élastique ;
the thermal term acts only on the right-hand side (equivalent thermal load
`f_th = ∫ Bᵀ D ε_th`) and on the real stress
`σ = D:(ε(u) − ε_th)`. Uncoupled: the mechanics does not feed back on the thermics.
Bricks composed by hand (no "all-in-one" operator):
`interp_to_gauss` (T nodale → points de Gauss), `thermal_strain` (EPTH),
`integrate_behavior` (σ = D:ε), `internal_forces` (BSIG), `solve`, puis
`deformation` and a field subtraction for ε_mech = ε(u) − ε_th.
`alpha` travels through the material field (`material_field`) as an
**optional** component of the elasticity, beside `E`/`nu`.
Two regimes on the same bar
------------------------------
- **bloquée** (deux bords en x encastrés) : `σ_xx = −E·α·ΔT`, `σ_yy = 0` ;
- **libre** (appuis simples) : dilatation `u = α·ΔT·(x, y)`, `σ ≈ 0`.
Lancement ::
maturin develop --features extension-module
python examples/thermoelastique_barre.py
"""
import pyrucast
E, NU, ALPHA = 210_000.0, 0.3, 1e-5
T_REF, DT = 20.0, 100.0
NX, NY, L, H = 4, 2, 4.0, 1.0
def _clamp(target, nodes, var):
imposed = pyrucast.mesh.poi1_from_nodes(nodes)
multiplier = pyrucast.mesh.barycenter(imposed)
return pyrucast.model.dirichlet(target, var, imposed, multiplier)
def _bar():
"""An NX×NY grid of QUA4 on [0,L]×[0,H]. Returns (coords, grid, fes, idx)."""
c = pyrucast.Coords(2)
hx, hy = L / NX, H / NY
def idx(i, j):
return j * (NX + 1) + i
grid = [c.add_node([i * hx, j * hy]) for j in range(NY + 1) for i in range(NX + 1)]
mesh = pyrucast.Mesh(c, "QUA4")
for j in range(NY):
for i in range(NX):
mesh.unit().add_cell(
[
grid[idx(i, j)],
grid[idx(i + 1, j)],
grid[idx(i + 1, j + 1)],
grid[idx(i, j + 1)],
]
)
return c, grid, pyrucast.FiniteElementSpace(mesh), idx
def _uniform_temperature(c, grid, fes, value):
"""A temperature field 'T' = value everywhere, carried at the Gauss points."""
t_mesh = pyrucast.Mesh(c, "POI1")
for node in grid:
t_mesh.unit().add_cell([node])
t_nodal = pyrucast.NodeField(t_mesh, ["T"])
for node in grid:
t_nodal[0].set_value(node, "T", value)
return pyrucast.element_field.interp_to_gauss(t_nodal, fes)
def _displacement(solution, c, grid):
"""Extracts a clean (u_x, u_y) field (without the Lagrange multipliers)."""
u_mesh = pyrucast.Mesh(c, "POI1")
for node in grid:
u_mesh.unit().add_cell([node])
u = pyrucast.NodeField(u_mesh, ["u_x", "u_y"])
for node in grid:
u[0].set_value(node, "u_x", solution.value(node, "u_x"))
u[0].set_value(node, "u_y", solution.value(node, "u_y"))
return u
def _solve_thermal(model, materials, fes, c, grid):
"""ε_th → charge thermique → u → σ = D:(ε(u) − ε_th)."""
eps_th = pyrucast.element_field.thermal_strain(
_uniform_temperature(c, grid, fes, T_REF + DT), materials, fes, T_REF
)
sig_th = pyrucast.element_field.integrate_behavior(model, eps_th, materials)
# A **load**, not a residual: the divergence of the prescribed tensor, hence
# the geometric operator. They are then renamed into dual rows — that is
# where, and only where, these numbers become forces.
f_th = (
pyrucast.node_field.divergence(sig_th, "sigma")
.rename_component("div_sigma_x", "f_x")
.rename_component("div_sigma_y", "f_y")
)
solution = pyrucast.solver.solve(pyrucast.matrix.stiffness(model, materials), f_th)
u = _displacement(solution, c, grid)
sigma = pyrucast.element_field.integrate_behavior(
model, pyrucast.element_field.deformation(u, fes) - eps_th, materials
)
return u, sigma
def main() -> None:
# ── Régime bloqué : σ_xx = −E·α·ΔT ──────────────────────────────────────
c, grid, fes, idx = _bar()
left = [grid[idx(0, j)] for j in range(NY + 1)]
right = [grid[idx(NX, j)] for j in range(NY + 1)]
bottom = [grid[idx(i, 0)] for i in range(NX + 1)]
model = pyrucast.model.elasticity(fes, "plane_stress")
model = (
model
| _clamp(model, left, "u_x")
| _clamp(model, right, "u_x")
| _clamp(model, bottom, "u_y")
)
materials = pyrucast.element_field.material_field(
model, [("E", E), ("nu", NU), ("alpha", ALPHA)]
)
_u, sigma = _solve_thermal(model, materials, fes, c, grid)
expected = -E * ALPHA * DT
sub = sigma[0]
sxx = sub.value(0, 0, "sigma_xx")
print(f"Bloquée : σ_xx = {sxx:12.4f} (attendu {expected:.4f} = −E·α·ΔT)")
assert abs(sxx - expected) < 1e-6 * abs(expected)
# ── Régime libre : dilatation u = α·ΔT·(x, y), σ ≈ 0 ────────────────────
c, grid, fes, idx = _bar()
left = [grid[idx(0, j)] for j in range(NY + 1)]
bottom = [grid[idx(i, 0)] for i in range(NX + 1)]
model = pyrucast.model.elasticity(fes, "plane_stress")
model = model | _clamp(model, left, "u_x") | _clamp(model, bottom, "u_y")
materials = pyrucast.element_field.material_field(
model, [("E", E), ("nu", NU), ("alpha", ALPHA)]
)
u, sigma = _solve_thermal(model, materials, fes, c, grid)
tip = grid[idx(NX, NY)]
ux, uy = u.value(tip, "u_x"), u.value(tip, "u_y")
print(
f"Libre : u(coin) = ({ux:.6e}, {uy:.6e}) (attendu ({ALPHA * DT * L:.6e}, {ALPHA * DT * H:.6e}))"
)
assert abs(ux - ALPHA * DT * L) < 1e-9 and abs(uy - ALPHA * DT * H) < 1e-9
assert abs(sigma[0].value(0, 0, "sigma_xx")) < 1e-6
print("\nOK: blocked bar → σ_xx = −E·α·ΔT; free bar → expansion without stress.")
if __name__ == "__main__":
main()