Plasticité parfaite (von Mises)
Élastoplasticité parfaite (sans écrouissage) en petites déformations,
critère de von Mises (J2), écoulement associé. Mêmes éléments et mêmes
degrés de liberté que l’élasticité linéaire : 2-D (TRI3 /
QUA4) ou 3-D (TET4 / HEX8).
Équations continues résolues
- équilibre :
∇·σ + b = 0; - partition :
ε = εᵉ + εᵖ(déformation élastique + plastique) ; - élasticité :
σ = D : εᵉ; - critère :
f(σ) = q − σ_y ≤ 0, avecq = √(3/2 · s:s)la contrainte équivalente de von Mises (s= déviateur deσ) ; - écoulement associé :
ε̇ᵖ = λ̇ · ∂f/∂σ, conditions de Kuhn–Tuckerλ̇ ≥ 0, f ≤ 0, λ̇ f = 0.
Sans écrouissage, σ_y est constant : la contrainte équivalente est
plafonnée à σ_y (plateau parfaitement plastique).
Forme discrétisée
Conformément au découpage du cœur (voir Comportement — la boucle de Newton est pilotée en Python, le cœur Rust ne fournit que les briques), cette physique expose :
- rigidité (
build_stiffness_blocks) : la rigidité élastiqueK = ∫ Bᵀ D B dΩ, opérateur d’itération simple (Newton modifié) ; - tangente cohérente (
KTAN,assemble.tangent) :K_t = ∫ Bᵀ D_alg Bavec le module algorithmiqueD_algdu retour radial J2 (dérivée exacte, condensation contrainte-plane incluse) — émis par le comportement, relu par l’assembleur, pour un Newton complet à convergence quadratique ; - comportement (
COMP,integrate_behavior) : le retour radial exact, point de Gauss par point de Gauss, qui produit aussiD_alg.
À chaque itération, avec la même matrice \( B \) qu’en élasticité, on résout la correction \( \delta u \) du système linéarisé
\[
K_t\,\delta u = F_{\text{ext}} - F_{\text{int}}, \qquad
F_{\text{int}} = \int_\Omega B^\top \sigma\, d\Omega \;(\text{BSIG}),
\]
où la contrainte \( \sigma \) au point de Gauss sort du retour radial. La tangente cohérente \( K_t = \int_\Omega B^\top D_{\text{alg}} B\, d\Omega \) utilise le module algorithmique \( D_{\text{alg}} = \partial\sigma/\partial\varepsilon \) (dérivée exacte de l’application de retour), garantissant la convergence quadratique — au lieu de la rigidité élastique \( D \) du Newton modifié.
Retour radial (algorithme)
À partir de la déformation totale ε et de l’état plastique précédent
(εᵖ, p) :
- prédiction élastique :
σ_trial = D : (ε − εᵖ),q = √(3/2 s_trial:s_trial); - si
f = q − σ_y ≤ 0→ pas élastique, état inchangé ; - sinon (plasticité parfaite) :
Δp = f / (3μ), le déviateur est ramenés = s_trial · σ_y / q,Δεᵖ = Δp · (3/2) s_trial / q, puisεᵖ ← εᵖ + Δεᵖ,p ← p + Δp.
Le calcul interne est mené en 3-D quel que soit le modèle. La
déformation plane impose ε_zz = ε_yz = ε_xz = 0 ; la contrainte plane
résout la condition σ_zz = 0 par une méthode de la sécante locale autour du
retour radial.
Variables et matériau
- primal :
u_x, u_y(, u_z)— dual :f_x, f_y(, f_z). - matériau :
E(Young),nu(Poisson),sigma_y(limite d’élasticité). - état de début de pas A (entrée
prev, montage incrémental) : contrainteσ(A), tenseur de déformation plastiqueeps_p_xx … eps_p_xy(toujours 6 composantes 3-D), déformation plastique cumuléep, et déformationε(A). C’est la sortie du pas précédent ;Noneau premier pas (A = configuration de référence, tout à zéro). - sortie du
COMP(état de B, =prevdu pas suivant) : contrainte (sigma_*dans l’ordre de Voigt du modèle) suivie de l’état mis à jour et de l’écho deε(B)full-3-D (plussigma_zzpour les modèles plans en 2-D, dont le dual de Voigt l’omet) pour queprevsoit complet. Le prédicteur élastique estσ_trial = σ(A) + C:(ε(B) − ε(A)). - modèles :
plane_stress,plane_strain,axisymmetric(2-D) etfull_3d(3-D).
Axisymétrie
Le modèle "axisymmetric" s’applique sur une géométrie de révolution
(Coords.axisymmetric()) : Voigt à quatre
composantes [εrr, εzz, εθθ, γrz], nommées eps_xx, eps_yy, eps_zz, eps_xy
avec zz = orthoradial (convention Cast3M). Le modèle et le repère doivent
s’accorder dans les deux sens, comme en élasticité.
L’état interne étant déjà stocké en 3-D complet et le retour radial se
faisant en 3-D, la spécialisation axisymétrique se réduit à une table d’indices
[rr, zz, θθ, rz] → [xx, yy, zz, xy]. Deux conséquences :
- la déformation orthoradiale
ε_θθ = u_r/rest mesurée (produite pardeformation), pas supposée :ε(B)est donc entièrement connue, sans la résolution hors-plan qu’exige la contrainte plane ; σ_zzfait partie du dual de Voigt et n’est donc pas ré-émis en écho, contrairement aux modèles plans.
La tangente cohérente axisymétrique est la restriction [rr, zz, θθ, rz] de la
tangente 3-D, validée par différences finies dans tests/tangent.rs.
Mise en donnée (Rust) : poutre console, boucle de Newton complète
examples/plasticite_poutre_console.rs déroule un Newton complet (et non
un seul pas) autour des briques ci-dessus, côté API Rust : une poutre encastrée
cisaillée au bout, chargée par incréments jusqu’à développer une zone plastique
à l’encastrement.
L’algorithmie de Newton vit entièrement dans l’exemple, pas dans pyrucast :
la bibliothèque ne fournit que les opérateurs ponctuels — stiffness (rigidité
élastique, opérateur d’itération), deformation (ε), integrate (COMP,
retour radial → σ + état), internal_forces (BSIG, ∫ Bᵀσ) et solve
(LU creux, factorisation en cache). L’exemple assemble lui-même le résidu
r = F_ext − F_int, résout δu = K⁻¹ r et porte l’état interne VAR0 → VAR1
d’un pas au suivant. C’est un Newton modifié : K élastique constant,
assemblé et factorisé une seule fois.
Étant en Rust pur (aucune dépendance à Python), il sert aussi de banc de
parallélisme — les boucles chaudes (assemblage, deformation, integrate,
internal_forces) sont réévaluées à chaque itération :
RAYON_NUM_THREADS=1 PYRUCAST_NX=200 PYRUCAST_NY=40 \
cargo run --release --example plasticite_poutre_console
Variables d’environnement : PYRUCAST_NX / PYRUCAST_NY (mailles),
PYRUCAST_NSTEPS (pas de charge), PYRUCAST_PMAX (charge finale).
Exemple Python
La boucle de Newton (assemblage du résidu, résolution, mise à jour de l’état) s’écrit en Python ; voici l’usage d’un pas de la brique d’intégration :
import pyrucast
model = pyrucast.model.plasticity_perfect(fes, "plane_stress")
materials = pyrucast.element_field.material_field(
model, [("E", 210_000.0), ("nu", 0.3), ("sigma_y", 250.0)]
)
# Strain ε(B) from the current displacement field (a geometric op).
strain = pyrucast.element_field.deformation(u, fes)
# Integration A→B: `prev` = the previous step's output (None at the first step).
state = pyrucast.element_field.integrate_behavior(
model, strain, materials, prev=prev_state
)
sigma_xx = state[0].value(0, 0, "sigma_xx")
p = state[0].value(0, 0, "p") # déformation plastique cumulée
Pour réinjecter l’état au pas suivant, il suffit de passer state comme
prev au prochain appel — la sortie porte déjà l’état complet de B (σ, VAR1,
ε(B)). Aucune fusion de champs n’est nécessaire. La boucle de Newton
complète (pas de charge, résidu, résolution, portage de l’état) est écrite
dans examples/plasticite_poutre_console.py — même architecture que la mise en
donnée Rust ci-dessus.