Élasticité orthotrope et anisotrope

En bref

Un composite tissé, un bois, une tôle laminée ou un monocristal n’ont pas la même raideur dans toutes les directions. C’est l’option symmetry de l’élasticité, indépendante de la cinématique :

Les variables ne changent pas : seule la matrice de Hooke \(D\) change. Les lois non linéaires (plasticité, fluage, endommagement de Mazars et Damage TC) gardent une élasticité isotrope ; la loi SiC/SiC a ses propres axes.

Théorie

La loi de comportement est \(\sigma_{ij} = C_{ijkl}\,\varepsilon_{kl}\). Les symétries de \(\sigma\) et \(\varepsilon\) et l’existence d’un potentiel élastique réduisent \(C\) à une matrice \(6\times6\) symétrique en notation de Voigt : 21 constantes dans le cas général.

Orthotropie. Dans ses axes propres, la souplesse \(S = C^{-1}\) découple les termes normaux des cisaillements :

\[ S = \begin{bmatrix} 1/E_1 & -\nu_{21}/E_2 & -\nu_{31}/E_3 & & & \\ -\nu_{12}/E_1 & 1/E_2 & -\nu_{32}/E_3 & & & \\ -\nu_{13}/E_1 & -\nu_{23}/E_2 & 1/E_3 & & & \\ & & & 1/G_{23} & & \\ & & & & 1/G_{13} & \\ & & & & & 1/G_{12} \end{bmatrix}, \qquad \frac{\nu_{ji}}{E_j} = \frac{\nu_{ij}}{E_i} . \]

Les relations de réciprocité ne laissent que neuf constantes. \(S\) doit être définie positive : un jeu de constantes qui la rend singulière est refusé, plutôt que d’assembler une raideur fausse.

Les axes matériau sont donnés par des vecteurs du champ matériau : V1X, V1Y en 2-D (le second axe est la normale dans le plan, le troisième la direction hors du plan, l’orthoradiale en axisymétrie) ; V1X…V1Z et V2X…V2Z en 3-D (le troisième est \(V_1 \times V_2\)). Ils sont orthonormalisés en interne : \(V_2\) n’a besoin d’être qu’à peu près perpendiculaire à \(V_1\). Un \(V_1\) nul, ou un \(V_2\) parallèle à \(V_1\), est refusé. Comme toute composante matériau, les axes peuvent varier d’une maille et d’un point de Gauss à l’autre : un composite bobiné, une pièce courbe.

En chaque point de Gauss, \(D\) est construit dans les axes matériau, tourné dans les axes globaux par le tenseur d’ordre 4,

\[ C'_{pqrs} = R_{pi}\,R_{qj}\,R_{rk}\,R_{sl}\,C_{ijkl}, \qquad R = \big[\,V_1\ \ V_2\ \ V_1\times V_2\,\big], \]

puis réduit à la cinématique : le bloc [xx, yy, xy] en déformations planes, sa condensation sur \(\varepsilon_{zz}\) en contraintes planes, le bloc [rr, zz, θθ, rz] en axisymétrie. En cisaillement de l’ingénieur, le passage Voigt ↔ tenseur ne porte aucun facteur, et la rotation d’ordre 4 n’a aucune convention à retenir.

Variables

Les neuf constantes orthotropes sont demandées même en 2-D : le tenseur est construit en 3-D avant d’être réduit. Les trois jeux de composantes sont disjoints, si bien qu’une zone isotrope et une zone orthotrope partagent un maillage sans rien consolider. L’entrée et la sortie du comportement sont celles de l’élasticité isotrope.

Exemple

Les constantes d’un matériau quatre fois plus raide selon son axe 1 :

gauche, droite, carre = _carre()
plaque = pyrucast.model.elasticity(
    pyrucast.FiniteElementSpace(carre), "plane_stress", symmetry="orthotropic"
)
print(plaque[0].material_components())
# ['E_1', 'E_2', 'E_3', 'nu_12', 'nu_13', 'nu_23', 'G_12', 'G_13', 'G_23', 'V1X', 'V1Y']
angle = 0.0  # l'axe raide, selon x
constantes = [
    ("E_1", 40e3),
    ("E_2", 10e3),
    ("E_3", 10e3),
    ("nu_12", 0.25),
    ("nu_13", 0.25),
    ("nu_23", 0.3),
    ("G_12", 5e3),
    ("G_13", 5e3),
    ("G_23", 4e3),
    ("V1X", math.cos(angle)),
    ("V1Y", math.sin(angle)),
]

Sur la plaque en traction de la page Élasticité, le fichier de tests vérifie les deux bornes : \(u_x = S/E_1\) quand l’axe raide suit la traction, \(S/E_2\) quand il lui est perpendiculaire.

Côté Rust

src/models/symmetry.rs construit les tenseurs orientés pour l’élasticité, la conduction et la diffusion ; MaterialSymmetry est l’énumération des trois symétries. Le test ci-dessous étire un carré orthotrope selon son premier axe matériau ; le même fichier nourrit une loi orthotrope de constantes isotropes, puis une loi anisotrope du tenseur isotrope, qui doivent se comporter comme l’isotrope quel que soit le repère, et refuse un repère dégénéré.

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;
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::symmetry::MaterialSymmetry;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::StrainMeasure;
use pyrucast::ops::model;
use pyrucast::ops::solver::solve;
use pyrucast::Result;

/// Traction on the right edge.
const S: f64 = 2.0;
/// `N×N` QUA4 grid on the unit square.
const N: usize = 2;

#[test]
fn orthotropic_square_stretches_along_its_first_material_axis() -> Result<()> {
    const E1: f64 = 200.0; // stiff direction — aligned with global x
    const E2: f64 = 50.0; // compliant transverse direction
    const NU12: f64 = 0.25;

    let (grid, fes, _coords) = unit_square()?;
    let model = clamped_model(&grid, &fes, MaterialSymmetry::Orthotropic)?;

    // The material axes travel through the material field like any other
    // coefficient: `V1` is the first orthotropy direction, here the global x.
    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[
            ("E_1", E1),
            ("E_2", E2),
            ("E_3", E2),
            ("nu_12", NU12),
            ("nu_13", NU12),
            ("nu_23", 0.25),
            ("G_12", 30.0),
            ("G_13", 30.0),
            ("G_23", 30.0),
            ("V1X", 1.0),
            ("V1Y", 0.0),
            ("phi_f_x", S),
        ],
    )?;

    let solution = solve_traction(&model, &materials)?;

    // σ_xx = S with σ_yy = σ_zz = 0 ⇒ the compliance row gives ε directly.
    let tol = 1e-10;
    let h = 1.0 / N as f64;
    for j in 0..=N {
        for i in 0..=N {
            let (x, y) = (i as f64 * h, j as f64 * h);
            let id = grid[j * (N + 1) + i].id();
            let ux = solution.value(id, "u_x")?;
            let uy = solution.value(id, "u_y")?;
            assert!((ux - S / E1 * x).abs() < tol, "u_x({x},{y}) = {ux}");
            assert!((uy + NU12 * S / E1 * y).abs() < tol, "u_y({x},{y}) = {uy}");
        }
    }
    Ok(())
}