Élasticité orthotrope et anisotrope

Introduction

L’élasticité linéaire suppose le matériau isotrope : deux constantes, aucune direction privilégiée. Beaucoup de matériaux n’obéissent pas à cette hypothèse — un composite tissé, un bois, une tôle laminée, un monocristal ont des raideurs différentes selon la direction.

pyrucast traite cela comme un axe à part entière, la symétrie matériau, porté par le même modèle elasticity et orthogonal à l’hypothèse cinématique (contraintes planes, déformations planes, axisymétrie, massif) :

C’est la convention de Cast3M, où ISOTROPE / ORTHOTROPE / ANISOTROPE qualifie le matériau d’une formulation plutôt que de nommer un modèle différent. Les degrés de liberté ne changent donc pas : déplacement u_x, u_y(, u_z) en primal, force nodale f_x, … en dual, exactement comme en isotrope. Seule la matrice de Hooke D change.

Équations continues résolues

Les mêmes qu’en élasticité linéaire — équilibre ∇·σ + f = 0 en petites déformations, avec ε = ½(∇u + ∇uᵀ). C’est la loi de comportement qui se généralise :

\[ \sigma_{ij} = C_{ijkl}\,\varepsilon_{kl} \]

où \( C \) est le tenseur d’élasticité d’ordre 4. Ses symétries — mineures \( C_{ijkl} = C_{jikl} = C_{ijlk} \), qui viennent de celles de \( \sigma \) et \( \varepsilon \), et majeure \( C_{ijkl} = C_{klij} \), qui vient de l’existence d’un potentiel élastique \( W = \tfrac12\,\varepsilon : C : \varepsilon \) — le réduisent d’un tenseur à 81 composantes à une matrice \( 6\times6 \) symétrique en notation de Voigt, soit 21 constantes dans le cas général.

L’orthotropie est le cas où le matériau possède trois plans de symétrie orthogonaux. Dans ses axes propres, la souplesse S = C⁻¹ se découple : les termes normaux ne sont couplés qu’entre eux, et chaque cisaillement est isolé.

\[ 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} \]

Les relations de réciprocité \( \nu_{ji}/E_j = \nu_{ij}/E_i \) rendent la matrice symétrique, d’où neuf constantes seulement : trois modules d’Young, trois coefficients de Poisson, trois modules de cisaillement. Le bloc normal et les trois cisaillements sont découplés, ce qui est la définition même de l’orthotropie : une traction selon un axe propre ne produit aucun cisaillement.

Attention — un jeu de constantes n’est pas physique par construction : S doit rester définie positive, ce qui impose ν_ij² < E_i/E_j. pyrucast le vérifie en inversant S et erronne si elle est singulière, plutôt que d’assembler une raideur non définie positive en silence.

Le repère d’orthotropie

Les constantes sont données dans les axes matériau, qui ne coïncident pas avec les axes globaux. Il faut donc dire où ils pointent.

pyrucast suit Cast3M et les décrit par des vecteurs, pas par des angles d’Euler (MATE 'DIRECTION' V1 V2). Ils voyagent dans le champ matériau comme n’importe quel autre coefficient :

Ils sont orthonormalisés en interne (Gram-Schmidt) : V2 n’a besoin d’être que grossièrement perpendiculaire à V1, c’est le plan qu’ils engendrent qui compte. Des vecteurs plutôt que des angles, parce qu’il n’y a aucune convention à retenir, aucun cas de blocage de cardan — et surtout parce que le repère varie alors naturellement d’une maille à l’autre (un composite bobiné, une pièce courbe), en passant par le canal matériau existant.

Un V1 nul, ou un V2 parallèle à V1, est un repère dégénéré : il est signalé, jamais complété arbitrairement.

Forme discrétisée

La chaîne est celle de l’élasticité — K_e = Σ_g Bᵀ D B |J| w, avec le même opérateur B. Ce qui change tient en trois étapes, faites une fois par maille :

1. construire D dans les axes matériau, où l’orthotropie est diagonale ;
2. le tourner vers les axes globaux ;
3. le réduire au modèle cinématique (bloc [xx, yy, xy] en déformations planes, sa condensation statique sur ε_zz en contraintes planes, le bloc [rr, zz, θθ, rz] en axisymétrie, le 6×6 complet en massif).

La rotation passe par le tenseur d’ordre 4, pas par une matrice de Bond 6×6 :

\[ 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], \]

R étant la rotation qui porte les axes matériau sur les axes globaux.

C’est un choix délibéré. En cisaillement ingénieur (γ = 2ε), le passage Voigt ↔ tenseur ne porte aucun facteur — C_ijkl = D[voigt(i,j)][voigt(k,l)] — et la rotation d’ordre 4 s’écrit sans la moindre convention à mémoriser. Le coût, quelques centaines de multiplications par maille, est négligeable devant l’assemblage, et il achète l’élimination de toute une famille d’erreurs d’indices et de facteurs 2. L’isotropie, elle, court-circuite ce chemin et garde ses formes fermées : les calculs isotropes existants conservent leurs nombres exacts.

Variables et matériau

Primales u_x, u_y(, u_z), duales f_x, f_y(, f_z) — inchangées.

Les trois contrats sont disjoints. Ce n’est pas un détail : l’assembleur résout la zone matériau d’une physique par l’ensemble des composantes qu’elle exige, si bien qu’une zone isotrope et une zone orthotrope peuvent partager un maillage sans consolidation explicite.

Les constantes anisotropes sont nommées d’après le triangle supérieur de la matrice de Voigt, dans l’ordre de ce dépôt [xx, yy, zz, yz, xz, xy] : C_11, C_12, …, C_16, C_22, …, C_66. Ainsi C_44 est la raideur en yz, C_66 celle en xy.

Même en 2-D, les neuf constantes orthotropes sont exigées : la raideur hors-plan intervient en déformations planes et en axisymétrie, et le tenseur complet est de toute façon construit avant d’être réduit.

Le comportement (COMP) est linéaire, σ = D·ε, et rend les mêmes composantes sigma_* que l’élasticité isotrope.

Mise en donnée (Rust, 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;
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::symmetry::MaterialSymmetry;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::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(())
}

Exemple Python

Le balayage du repère matériau, où l’on voit l’effet propre à l’orthotropie :

"""Orthotropic elasticity — a plate pulled off its material axes.

Problème
--------
A unit square in plane stress, pulled uniformly (traction ``S``) on its right
edge, with roller supports on the left and bottom edges. The material is
**orthotropic**: stiff in one direction, compliant in the other.

What the example shows is the effect of the **orthotropy frame**. It is given
by vectors, as in Cast3M (``MATE 'DIRECTION' V1 V2``): the components ``V1X``,
``V1Y`` travel in the material field just like the moduli. The first material
axis's angle is swept from 0° to 90°.

Two cases have an analytical solution, and they are the sweep's bounds:

  * **0°** — the stiff axis is aligned with the traction ::

        u_x(1, y) = S / E_1

  * **90°** — it is the compliant axis that works ::

        u_x(1, y) = S / E_2

In between, the plate **shears**: off its axes, an orthotropic material couples
traction and distortion (the ``D_16`` term of the rotated tensor is no longer
zero), and the right edge does not stay straight. That is precisely what
anisotropy brings, and what an isotropic computation cannot produce.

Lancement
---------
Once the extension is built in the venv ::

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

import math

import pyrucast

# ── Problem data ────────────────────────────────────────────────────────────
E1 = 200.0  # modulus in the stiff direction (material axis 1)
E2 = 50.0  # module transverse
NU12 = 0.25
G12 = 30.0
S = 2.0  # traction on the right edge
N = 4  # grille N×N de QUA4


def maillage():
    """The unit square's QUA4 grid, its nodes and its FE space."""
    h = 1.0 / N
    c = pyrucast.Coords(2)
    grid = [[c.add_node([i * h, j * h]) for i in range(N + 1)] for j in range(N + 1)]
    mesh = pyrucast.Mesh(c, "QUA4")
    for j in range(N):
        for i in range(N):
            mesh.unit().add_cell(
                [grid[j][i], grid[j][i + 1], grid[j + 1][i + 1], grid[j + 1][i]]
            )
    return c, grid, pyrucast.FiniteElementSpace(mesh)


def rouleau(target, c, noeuds, variable):
    """A roller support ``variable = 0`` on the given nodes."""
    imposed = pyrucast.Mesh(c, "POI1")
    for n in noeuds:
        imposed.unit().add_cell([n])
    multiplier = pyrucast.mesh.barycenter(imposed)
    return pyrucast.model.dirichlet(target, variable, imposed, multiplier)


def resoudre(angle_deg):
    """The displacement of corner (1, 0) for a material axis at ``angle_deg``."""
    c, grid, fes = maillage()

    # Orthotropic elasticity + both supports.
    model = pyrucast.model.elasticity(fes, "plane_stress", symmetry="orthotropic")
    model = model | rouleau(model, c, [grid[j][0] for j in range(N + 1)], "u_x")
    model = model | rouleau(model, c, [grid[0][i] for i in range(N + 1)], "u_y")

    # The material frame is material data like any other.
    a = math.radians(angle_deg)
    # Traction S on the right edge, as consistent nodal loads: a term of the
    # model, whose density lives in the material.
    bord = pyrucast.Mesh(c, "SEG2")
    for j in range(N):
        bord.unit().add_cell([grid[j][N], grid[j + 1][N]])
    bord_fes = pyrucast.FiniteElementSpace(bord)
    model = model | pyrucast.model.flux(bord_fes, model, "f_x")

    materials = pyrucast.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", G12),
            ("G_13", G12),
            ("G_23", G12),
            ("V1X", math.cos(a)),
            ("V1Y", math.sin(a)),
            ("phi_f_x", S),
        ],
    )
    rhs = pyrucast.node_field.external_forces(model, materials)

    solution = pyrucast.solver.solve(pyrucast.matrix.stiffness(model, materials), rhs)
    coin = grid[0][N]  # (1, 0)
    haut = grid[N][N]  # (1, 1)
    return (
        solution.value(coin, "u_x"),
        solution.value(haut, "u_x") - solution.value(coin, "u_x"),
    )


def main() -> None:
    print("Orthotropic elasticity — sweeping the material frame")
    print(f"  E_1 = {E1}, E_2 = {E2}, nu_12 = {NU12}, G_12 = {G12}, traction S = {S}")
    print()
    print("  angle    u_x(1,0)    u_x gap on the right edge")
    print("  " + "-" * 46)
    for angle in (0.0, 22.5, 45.0, 67.5, 90.0):
        ux, distorsion = resoudre(angle)
        print(f"  {angle:5.1f}°  {ux:10.6f}  {distorsion:+14.6f}")

    # Both bounds are analytical: the stiff axis, then the compliant one.
    ux0, _ = resoudre(0.0)
    ux90, _ = resoudre(90.0)
    print()
    print(f"  0°  : {ux0:.6f}  (attendu S/E_1 = {S / E1:.6f})")
    print(f"  90° : {ux90:.6f}  (attendu S/E_2 = {S / E2:.6f})")
    assert abs(ux0 - S / E1) < 1e-10
    assert abs(ux90 - S / E2) < 1e-10

    # Off axis, the traction induces shear — the right edge warps.
    _, distorsion45 = resoudre(45.0)
    assert abs(distorsion45) > 1e-4, "un orthotrope hors axes doit cisailler"
    print()
    print(f"  At 45°, the right edge warps by {distorsion45:+.6f}:")
    print("  this is the traction/shear coupling of off-axis orthotropy.")


if __name__ == "__main__":
    main()

Il produit :

  angle    u_x(1,0)    écart u_x sur le bord droit
  ----------------------------------------------
    0.0°    0.010000       -0.000000
   22.5°    0.017888       -0.004075
   45.0°    0.031195       -0.008035
   67.5°    0.039603       -0.005663
   90.0°    0.040000       -0.000000

Les deux bornes sont analytiques — S/E₁ quand l’axe rigide est aligné sur la traction, S/E₂ quand c’est l’axe souple. Entre les deux, le bord droit se gauchit : hors de ses axes, un matériau orthotrope couple traction et cisaillement (le terme D₁₆ du tenseur tourné n’est plus nul). C’est exactement ce qu’un calcul isotrope ne peut pas produire, et le signe le plus visible que la rotation fait son travail.

Compléments

Ce qui vaut comme vérification. Deux dégénérescences encadrent l’implémentation, et sont testées de bout en bout :

• une loi orthotrope nourrie de constantes isotropes doit se comporter comme l’isotrope, quel que soit son repère — c’est le contrôle le plus sévère de la rotation, puisque toute erreur d’indice brise l’invariance ;
• une loi anisotrope nourrie du tenseur isotrope doit faire de même, ce qui fixe l’ordre de lecture des 21 constantes : une permutation placerait les modules de cisaillement dans les mauvaises cases de Voigt.

Axisymétrie. L’orthotropie s’y combine sans rien de particulier. En 2-D le troisième axe matériau est la direction hors-plan, c’est-à-dire l’orthoradiale θ — ce qui est le comportement voulu pour un tube bobiné, dont la direction de fibre est justement circonférentielle.

Ce que cela ne couvre pas. La symétrie matériau porte sur l’élasticité. Les lois non linéaires (plasticité, endommagement) restent bâties sur une élasticité isotrope ; l’endommagement orthotrope est le sujet d’un modèle propre, pas d’un axe de symétrie.

La même mécanique sert la conduction thermique orientée et la diffusion, avec un tenseur d’ordre 2 au lieu de 4.