Plasticité de von Mises

En bref

Un métal qui se déforme de façon irréversible au-delà d’une limite d’élasticité. Deux lois, au critère de von Mises :

Théorie

L’élastoplasticité en petites déformations partage la déformation en une part élastique et une part plastique :

\[ \varepsilon = \varepsilon^e + \varepsilon^p, \qquad \sigma = D : (\varepsilon - \varepsilon^p), \qquad f(\sigma, p) \le 0, \qquad \dot\varepsilon^p = \dot\lambda\,\frac{\partial f}{\partial\sigma}, \]

avec les conditions de Kuhn-Tucker \(\dot\lambda \ge 0\), \(f \le 0\), \(\dot\lambda\,f = 0\). Le critère de von Mises porte sur la contrainte équivalente \(q = \sqrt{\tfrac32\,s : s}\), où \(s\) est le déviateur de \(\sigma\) :

\[ f(\sigma, p) = q - \big(\sigma_y + H\,p\big), \qquad \dot p = \dot\lambda , \]

\(p\) étant la déformation plastique cumulée et \(H\) le module d’écrouissage, nul pour la loi parfaite.

Le retour radial. Un pas part de l’état \(A\) (\(\varepsilon^p_A\), \(p_A\)) et de la déformation totale \(\varepsilon_B\) :

1. prédiction élastique, \(\sigma^{tr} = D : (\varepsilon_B - \varepsilon^p_A)\), et sa contrainte équivalente \(q^{tr}\) ;
2. si \(q^{tr} \le \sigma_y + H\,p_A\), le pas est élastique ;
3. sinon, la normale étant colinéaire au déviateur, la consistance est affine et se résout sans itérer :

\[ \Delta p = \frac{q^{tr} - (\sigma_y + H\,p_A)}{3\mu + H}, \qquad s_B = s^{tr}\Big(1 - \frac{3\mu\,\Delta p}{q^{tr}}\Big), \qquad \varepsilon^p_B = \varepsilon^p_A + \frac{3\,\Delta p}{2\,q^{tr}}\,s^{tr}, \]

\(\mu\) étant le module de cisaillement ; la pression ne change pas.

Le calcul est mené en 3-D quelle que soit la cinématique. En déformations planes, les composantes hors du plan de la déformation sont nulles ; en contraintes planes, la condition \(\sigma_{zz} = 0\) est résolue par une itération locale autour du retour radial ; en axisymétrie, la déformation orthoradiale est mesurée.

La tangente cohérente est la dérivée exacte du retour,

\[ D_\text{alg} = K\,I\otimes I + 2\mu\,\theta\,\mathbb I_\text{dev} - 2\mu\,\bar\theta\;\hat n\otimes\hat n, \qquad \theta = \frac{\sigma_y + H\,p_B}{q^{tr}}, \qquad \bar\theta = \frac{3\mu}{3\mu + H} - (1 - \theta), \]

avec \(\hat n = s^{tr}/\lVert s^{tr}\rVert\) et \(K\) le module de compressibilité. Le facteur \(\theta\) rend compte du pas fini : c’est ce qui donne à Newton sa convergence quadratique.

La raideur reste la raideur élastique \(\int B^{\mathsf T} D\,B\) : un opérateur constant, celui d’un Newton modifié. La tangente \(\int B^{\mathsf T} D_\text{alg} B\) est celle d’un Newton complet.

Variables

Selon la cinématique :

En axisymétrie, eps_zz et sigma_zz sont les composantes orthoradiales.

La sortie du comportement est la contrainte, suivie de l’état :

La déformation plastique et la déformation totale ont toujours six composantes, quelle que soit la cinématique. Au premier pas, prev est l’état de repos, element_field.initial_state(modele, materiaux), où tout est nul. Ensuite, la sortie d’un pas se passe telle quelle comme prev au pas suivant.

Ce qu’elle produit

Exemple

Une maille étirée à \(\varepsilon_{xx} = 2\cdot10^{-3}\), au-delà de \(\sigma_y/E\). La contrainte équivalente reste à \(\sigma_y\) :

fes, u = _un_qua4(2e-3)  # une maille, ε_xx = 2·10⁻³ : au-delà de σ_y/E
modele = pyrucast.model.plasticity_perfect(fes, "plane_stress")
materiaux = pyrucast.element_field.material_field(
    modele, [("E", 210e3), ("nu", 0.3), ("sigma_y", 250.0)]
)

eps = pyrucast.element_field.strain(modele, u, materiaux)
repos = pyrucast.element_field.initial_state(
    modele, materiaux
)  # premier pas : l'état de repos
etat = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, repos)
s = {nom: etat[0].value(0, 0, nom) for nom in ("sigma_xx", "sigma_yy", "sigma_xy")}
q = math.sqrt(
    s["sigma_xx"] ** 2
    + s["sigma_yy"] ** 2
    - s["sigma_xx"] * s["sigma_yy"]
    + 3 * s["sigma_xy"] ** 2
)
print(round(q, 3))  # 250.0 : la contrainte équivalente est plafonnée à σ_y
print(etat[0].value(0, 0, "p") > 0)  # True : de la déformation plastique

# Le pas suivant repart de cet état.
etat = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, prev=etat)
Kt = pyrucast.matrix.tangent(
    modele, materiaux, eps, prev=etat
)  # le module algorithmique

Avec écrouissage, elle vaut \(\sigma_y + H\,p\) :

modele = pyrucast.model.plasticity_isotropic(fes, "plane_stress")
materiaux = pyrucast.element_field.material_field(
    modele, [("E", 210e3), ("nu", 0.3), ("sigma_y", 250.0), ("H", 1000.0)]
)
eps = pyrucast.element_field.strain(modele, u, materiaux)
etat = pyrucast.element_field.integrate_behavior(
    modele, eps, materiaux, pyrucast.element_field.initial_state(modele, materiaux)
)
s = {nom: etat[0].value(0, 0, nom) for nom in ("sigma_xx", "sigma_yy", "sigma_xy")}
q = math.sqrt(
    s["sigma_xx"] ** 2
    + s["sigma_yy"] ** 2
    - s["sigma_xx"] * s["sigma_yy"]
    + 3 * s["sigma_xy"] ** 2
)
print(round(q - (250.0 + 1000.0 * etat[0].value(0, 0, "p")), 3))  # 0.0 : q = σ_y + H·p

La boucle de Newton complète (pas de charge, résidu \(F_\text{ext} - F_\text{int}\), résolution, report de l’état) s’écrit en Python avec ces opérateurs ; examples/plasticite_poutre_console.py en est un exemple, une console cisaillée jusqu’à plastifier à l’encastrement.

Côté Rust

Plasticity (src/models/plasticity.rs) porte la loi d’écoulement en attribut : une même physique sert les dix lois de cette page, de la page Drucker-Prager et Ottosen, de la page Fluage et la loi de Gurson. Les deux lois de von Mises vivent dans src/models/plasticity/von_mises.rs ; H = 0 redonne la loi parfaite par le même chemin de code. examples/plasticite_poutre_console.rs est la même console en Rust, et sert de banc de parallélisme.

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::element_field::ElementField;
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::plasticity::law::PlasticLaw;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::{
    behavior::{initial_state, integrate},
    deformation, material_field, StrainMeasure,
};
use pyrucast::ops::matrix::tangent;
use pyrucast::ops::model;
use pyrucast::ops::node_field::internal_forces;
use pyrucast::Result;

const AXES: [&str; 3] = ["x", "y", "z"];

/// von Mises with linear isotropic hardening.
const ISOTROPIC: &[(&str, f64)] = &[
    ("E", 70_000.0),
    ("nu", 0.3),
    ("sigma_y", 200.0),
    ("H", 5_000.0),
];
/// A frictional, mildly dilatant material — `ψ < α`, so the flow is
/// non-associated.
const DRUCKER: &[(&str, f64)] = &[
    ("E", 20_000.0),
    ("nu", 0.2),
    ("friction", 0.3),
    ("k", 30.0),
    ("psi", 0.1),
];
/// Ottosen's classic concrete set, for a tensile/compressive strength ratio of
/// about 0.1.
const OTTOSEN: &[(&str, f64)] = &[
    ("E", 30_000.0),
    ("nu", 0.2),
    ("a", 1.2759),
    ("b", 3.1962),
    ("k_1", 11.7365),
    ("k_2", 0.9801),
    ("sigma_c", 30.0),
];

#[test]
fn isotropic_hardening_satisfies_its_consistency_condition() -> Result<()> {
    let cube = Cube::new(PlasticLaw::Isotropic, ISOTROPIC)?;
    // Well past yield, so the step is definitely plastic.
    let s = cube.stress(&uniaxial(0.02))?;
    let (q, p) = (von_mises(&s.sigma), s.p);
    assert!(p > 0.0, "the step must be plastic (p = {p})");
    let expected = 200.0 + 5_000.0 * p;
    assert!(
        (q - expected).abs() < 1e-6 * expected,
        "q = {q}, expected σ_y + H·p = {expected}"
    );
    Ok(())
}