Dilatation, matériau homogène

La chape de la page Élasticité, encastrée et chargée par la force surfacique verticale, est maintenant chaude : 250 °C à l’alésage, une vingtaine de degrés près de l’encastrement. C’est la température du dernier instant du calcul de la page Thermique non linéaire, relue dans le fichier qu’il a écrit. La pièce se dilate, l’encastrement l’en empêche en partie, et des contraintes naissent sans aucun effort de plus.

Le script complet est tutoriel/mecanique_dilatation.py ; tous les extraits ci-dessous en sont issus directement, dans l’ordre du fichier. Il relit le fichier thermique.pyr qu’écrit tutoriel/thermique_non_lineaire.py : ce script-là doit avoir tourné avant.

L’équation résolue

Au système de la page Élasticité s’ajoute le terme de dilatation de la page Mécanique linéaire :

\[ [K]\{U\} = \{F_s\} + \{F_{\text{th}}\} + \{F_{\text{reac}}\}, \qquad \{F_{\text{th}}\} = \int_V [B]^T [C]\,\varepsilon_{\text{th}} \, dV, \qquad \varepsilon_{\text{th}} = \alpha\,(T - T_{\text{ref}})\,I \]

\( [C]\,\varepsilon_{\text{th}} \) est la contrainte qu’il faudrait pour empêcher toute dilatation ; \( \{F_{\text{th}}\} \) en sont les forces nodales. Le système résolu est le même que pour l’élasticité, avec ce second membre :

\[ \begin{bmatrix} K & A^T \\ A & 0 \end{bmatrix} \begin{Bmatrix} U \\ \lambda \end{Bmatrix} = \begin{Bmatrix} F_s + F_{\text{th}} \\ 0 \end{Bmatrix} \]

Seule la part de la déformation que la dilatation n’explique pas crée une contrainte :

\[ \sigma = [C] \left( [B]\{U\} - \varepsilon_{\text{th}} \right) \]

La dilatation ne pousse pas la pièce dans son ensemble. Les forces de \( \{F_{\text{th}}\} \) se compensent : leur somme est nulle, quelle que soit la température. Une pièce libre se dilaterait sans bouger dans son ensemble ; l’encastrement, lui, la retient, et ses réactions changent nœud par nœud. Leur somme, elle, reste le poids.

Chaque terme correspond à un objet du script :

Les données du calcul

E, NU = 200e9, 0.3  # Pa, — : un acier / a steel
ALPHA = 1e-5  # 1/K : coefficient de dilatation / expansion coefficient
T_REF = 0.0  # °C : la pièce n'est pas dilatée à 0 °C / not expanded at 0 °C
MASSE, G = 2500.0, 9.81  # kg, m/s²
POIDS = MASSE * G  # N
PRESSION = POIDS / (math.pi * HOLE_RADIUS * THICKNESS)  # Pa
AMPLIFICATION = 50.0

# FR — Le fichier écrit par tutoriel/thermique_non_lineaire.py.
# EN — The file tutoriel/thermique_non_lineaire.py writes.
FICHIER = os.path.join(tempfile.gettempdir(), "thermique.pyr")

TOL = 1e-9  # m

Les données de la page Élasticité, plus le coefficient de dilatation d’un acier, \( \alpha = 10^{-5} \) /K, et la température \( T_{\text{ref}} = 0 \) °C à laquelle la pièce n’est pas dilatée. FICHIER est le chemin où tutoriel/thermique_non_lineaire.py a écrit son calcul, dans le dossier temporaire du système.

La température relue

    # FR — Le volume et l'évolution des températures, relus sous les clefs
    #      où ils ont été écrits ; la température du dernier instant, en °C.
    # EN — The volume and the temperature evolution, read back under the
    #      keys they were written under; the last instant's temperature, in °C.
    if not os.path.exists(FICHIER):
        raise SystemExit(
            f"{FICHIER} manque : lancer d'abord tutoriel/thermique_non_lineaire.py"
        )
    relu = pc.load(FICHIER)
    volume = relu["volume"]
    evolution = relu["température (°C)"]
    temperature = evolution.frames()[-1]
    peau = pc.mesh.consolidate(pc.mesh.skin(volume))
    volume.plot(
        view=VUE,
        title="Température relue, dernier instant (°C)",
        field=temperature,
        component="T",
        cmap="viridis",
        save=os.path.join(OUT, "mecanique-dilatation-temperature.svg") if OUT else None,
    )

pc.load rend un dictionnaire, avec les clefs sous lesquelles pc.save a écrit : le volume, et l’évolution des 101 champs de température. frames() en donne les champs dans l’ordre du temps, et le dernier est celui de 15 000 s. Le volume relu est celui du calcul thermique, avec ses propres nœuds : tout ce qui suit se construit sur lui, et la température s’y lit sans rien convertir.

Une thermique et une mécanique se chaînent par un fichier. Le calcul thermique est fait une fois ; la mécanique le relit autant de fois qu’il le faut, sans le refaire.

La température relue va de 21,67 à 250,00 °C : l’alésage à 250 °C, le coin bas de l’encastrement le plus froid, refroidi à la fois par le flux qui sort de la face gauche et par le film convectif du dessous.

Les régions et le modèle

    # FR — Les régions de la page « Élasticité », retrouvées sur le volume
    #      relu : l'encastrement et la moitié basse de l'alésage.
    # EN — The "Élasticité" page's regions, found again on the volume read
    #      back: the clamp and the lower half of the bore.
    encastrement = pc.mesh.points_on_plane(volume, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0])
    bas_axe = [LENGTH, -THICKNESS, HEIGHT / 2.0]
    haut_axe = [LENGTH, 2.0 * THICKNESS, HEIGHT / 2.0]
    noeuds_alesage = pc.mesh.points_on_cylinder(peau, bas_axe, haut_axe, HOLE_RADIUS)
    alesage = pc.mesh.elements_on(peau, noeuds_alesage, strict=True)
    z_alesage = pc.node_field.positions(alesage, ["Z"])
    noeuds_bas = pc.mesh.select(z_alesage, le=HEIGHT / 2.0 + TOL)
    demi_alesage = pc.mesh.consolidate(
        pc.mesh.elements_on(alesage, noeuds_bas, strict=True)
    )

Les régions de la page Élasticité, retrouvées sur le volume relu par les mêmes sélections géométriques.

    # FR — Le modèle de la page « Élasticité » ; le matériau porte en plus
    #      « alpha », que la dilatation lira.
    # EN — The "Élasticité" page's model; the material also carries "alpha",
    #      which the expansion will read.
    fes = pc.FiniteElementSpace(volume)
    elasticite = pc.model.elasticity(fes, "full_3d")
    blocage_x = pc.model.dirichlet(encastrement, elasticite, "u_x")
    blocage_y = pc.model.dirichlet(encastrement, elasticite, "u_y")
    blocage_z = pc.model.dirichlet(encastrement, elasticite, "u_z")
    modele = elasticite | blocage_x | blocage_y | blocage_z
    materiaux = pc.element_field.material_field(
        modele, [("E", E), ("nu", NU), ("alpha", ALPHA)]
    )
    K = pc.matrix.stiffness(modele, materiaux)

    demi_fes = pc.FiniteElementSpace(demi_alesage)
    verticale = pc.model.flux(demi_fes, elasticite, "f_z")
    densite_verticale = pc.element_field.material_field(
        verticale, [("phi_f_z", -PRESSION)]
    )
    force_surfacique = pc.node_field.external_forces(verticale, densite_verticale)

Le modèle, l’encastrement et la force surfacique sont ceux de la page Élasticité. Le champ matériau porte en plus alpha, que la dilatation va lire. L’élasticité, elle, ne lit que E et nu : alpha ne change pas \( [K] \).

La dilatation

    # FR — La température aux points de Gauss, puis la dilatation
    #      ε_th = α (T − T_ref), la contrainte qu'elle donnerait si rien ne la
    #      gênait, [C] ε_th, et ses forces nodales ∫ [B]ᵀ [C] ε_th dV.
    # EN — The temperature at the Gauss points, then the expansion
    #      ε_th = α (T − T_ref), the stress it would give were nothing in its
    #      way, [C] ε_th, and its nodal forces ∫ [B]ᵀ [C] ε_th dV.
    t_gauss = pc.element_field.interp_to_gauss(temperature, fes)
    eps_th = pc.element_field.thermal_strain(t_gauss, materiaux, fes, T_REF)
    sigma_th = pc.element_field.elastic_stress(modele, eps_th, materiaux)
    force_th = (
        pc.node_field.divergence(sigma_th, "sigma")
        .rename_component("div_sigma_x", "f_x")
        .rename_component("div_sigma_y", "f_y")
        .rename_component("div_sigma_z", "f_z")
    )

Quatre lignes, une par terme du tableau ci-dessus. La température vient d’un calcul aux nœuds ; interp_to_gauss la porte aux points de Gauss, là où thermal_strain calcule \( \alpha\,(T - T_{\text{ref}}) \). elastic_stress en fait \( [C]\,\varepsilon_{\text{th}} \), et divergence intègre \( \int [B]^T \sigma \, dV \) sur toute la pièce.

divergence nomme ses composantes d’après celles du champ qu’elle reçoit : div_sigma_x, div_sigma_y, div_sigma_z. Le second membre de l’élasticité attend f_x, f_y, f_z : rename_component les renomme.

Le second membre

    # FR — Les deux chargements partagent les nœuds de l'alésage : la force
    #      surfacique est ramenée sur le support de la dilatation, puis `+`
    #      somme nœud à nœud.
    # EN — Both loads share the bore's nodes: the surface force is brought
    #      onto the expansion's support, then `+` sums node by node.
    second_membre = force_th + pc.node_field.restrict_like(force_surfacique, force_th)

Les deux chargements touchent les mêmes nœuds, ceux de la moitié basse de l’alésage. Comme à la page Régime stationnaire, leurs contributions doivent s’y additionner, ce que | ne fait pas. restrict_like ramène la force surfacique sur le support de la dilatation, tous les nœuds de la pièce ; les deux champs ont alors le même support, et + les somme nœud par nœud.

Résolution

    u = pc.solver.solve(K, second_membre)
    r_x = -pc.node_field.restrict(u, blocage_x.multiplier_mesh()).sum("lambda")
    r_y = -pc.node_field.restrict(u, blocage_y.multiplier_mesh()).sum("lambda")
    r_z = -pc.node_field.restrict(u, blocage_z.multiplier_mesh()).sum("lambda")

Les réactions se lisent comme à la page Élasticité. Le script affiche une somme des forces de la dilatation de l’ordre de \( 10^{-9} \) N, zéro à l’arrondi près, et des réactions de 0, 0 et 24 525,0 N : exactement le poids, comme sans dilatation.

La déformée

    coords = volume.coords()
    deformee = coords.add_config("déformée ×50")
    pc.coords.displace_from(u * AMPLIFICATION, 0, deformee)
    with coords.in_config(deformee):
        volume.plot(
            view=VUE,
            title="Dilatation et force surfacique — déformée ×50, u_x (mm)",
            field=u * 1000.0,
            component="u_x",
            cmap="viridis",
            save=os.path.join(OUT, "mecanique-dilatation-deformee.svg")
            if OUT
            else None,
        )

La pièce s’allonge : le bout libre avance de 0,760 mm selon \( x \). Chaque tranche s’allonge selon sa température, et les allongements s’additionnent de l’encastrement jusqu’au bout libre. Elle descend aussi davantage : le point le plus bas est à −1,280 mm, contre −0,978 mm sous la seule force surfacique. Le dessous, refroidi par le film convectif, se dilate moins que le dessus : la pièce se cintre vers le bas, dans le même sens que la masse.

Les contraintes

    # FR — Seule la part de la déformation que la dilatation n'explique pas
    #      crée une contrainte : σ = [C] (ε − ε_th).
    # EN — Only the part of the strain the expansion does not account for
    #      creates a stress: σ = [C] (ε − ε_th).
    deformation = pc.element_field.strain(modele, u, materiaux)
    contrainte = pc.element_field.elastic_stress(
        modele, deformation - eps_th, materiaux
    )
    von_mises_mpa = pc.element_field.von_mises(modele, contrainte) * 1e-6
    volume.plot(
        view=VUE,
        title="Dilatation et force surfacique — von Mises (MPa)",
        field=von_mises_mpa,
        component="von_mises",
        cmap="viridis",
        save=os.path.join(OUT, "mecanique-dilatation-von-mises.svg") if OUT else None,
    )

Aux points de Gauss, von Mises va de 3,92 à 379,19 MPa, contre 226,69 MPa au plus sans dilatation. Le maximum est au coin bas de l’encastrement : la face \( x = 0 \) ne peut pas se dilater, et le bas y est en plus comprimé par la flexion.

Un calcul élastique peut dépasser la limite du matériau. 379 MPa, c’est plus que les 250 MPa à partir desquels l’acier de la page Plasticité se déforme de façon permanente. Le calcul élastique ne le sait pas : il prolonge la loi de Hooke aussi loin qu’il le faut. Ce dépassement dit seulement que la pièce, en réalité, plastifierait.

Script complet

"""Tutoriel débutant — 9.7.2 Dilatation, matériau homogène. / Beginner
tutorial — 9.7.2 Thermal expansion, homogeneous material.

FR — Reprend la chape encastrée de la page « Élasticité », chargée par la
force surfacique verticale, et y ajoute la **dilatation thermique**. La
température est celle du dernier instant du calcul avec rayonnement, relue
dans le fichier qu'écrit `tutoriel/thermique_non_lineaire.py` : ce script-là
doit avoir tourné avant. La pièce se dilate de `α (T − T_ref)`, l'encastrement
l'en empêche en partie, et des contraintes naissent. Le détail pas à pas est
dans le livre, page « Dilatation, matériau homogène ».

EN — Takes the clamped lug of the "Élasticité" page back, loaded by the
vertical surface force, and adds **thermal expansion**. The temperature is
the last instant of the analysis with radiation, read back from the file
`tutoriel/thermique_non_lineaire.py` writes: that script must have run
first. The part expands by `α (T − T_ref)`, the clamp partly prevents it,
and stresses arise. The step-by-step walkthrough lives in the book's
"Dilatation, matériau homogène" page.

Lancement / Running ::

    pip install pyrucast
    python tutoriel/thermique_non_lineaire.py   # écrit / writes thermique.pyr
    python tutoriel/mecanique_dilatation.py

    # Figures du livre / book figures (book/src/tutoriel/img/) :
    # PYRUCAST_TUTORIEL_IMG_DIR=book/src/tutoriel/img python tutoriel/mecanique_dilatation.py
"""

import math
import os
import tempfile

import pyrucast as pc
from maillage import HEIGHT, HOLE_RADIUS, LENGTH, OUT, THICKNESS

# FR — La vue des figures de tutoriel/maillage.py : azimut, élévation, échelle.
# EN — The view of tutoriel/maillage.py's figures: azimuth, elevation, scale.
VUE = (-45, 25, 1.0)

# ── Données physiques / Physical data ──────────────────────────────────────
E, NU = 200e9, 0.3  # Pa, — : un acier / a steel
ALPHA = 1e-5  # 1/K : coefficient de dilatation / expansion coefficient
T_REF = 0.0  # °C : la pièce n'est pas dilatée à 0 °C / not expanded at 0 °C
MASSE, G = 2500.0, 9.81  # kg, m/s²
POIDS = MASSE * G  # N
PRESSION = POIDS / (math.pi * HOLE_RADIUS * THICKNESS)  # Pa
AMPLIFICATION = 50.0

# FR — Le fichier écrit par tutoriel/thermique_non_lineaire.py.
# EN — The file tutoriel/thermique_non_lineaire.py writes.
FICHIER = os.path.join(tempfile.gettempdir(), "thermique.pyr")

TOL = 1e-9  # m


def main() -> None:
    # FR — Le volume et l'évolution des températures, relus sous les clefs
    #      où ils ont été écrits ; la température du dernier instant, en °C.
    # EN — The volume and the temperature evolution, read back under the
    #      keys they were written under; the last instant's temperature, in °C.
    if not os.path.exists(FICHIER):
        raise SystemExit(
            f"{FICHIER} manque : lancer d'abord tutoriel/thermique_non_lineaire.py"
        )
    relu = pc.load(FICHIER)
    volume = relu["volume"]
    evolution = relu["température (°C)"]
    temperature = evolution.frames()[-1]
    peau = pc.mesh.consolidate(pc.mesh.skin(volume))
    volume.plot(
        view=VUE,
        title="Température relue, dernier instant (°C)",
        field=temperature,
        component="T",
        cmap="viridis",
        save=os.path.join(OUT, "mecanique-dilatation-temperature.svg") if OUT else None,
    )

    # FR — Les régions de la page « Élasticité », retrouvées sur le volume
    #      relu : l'encastrement et la moitié basse de l'alésage.
    # EN — The "Élasticité" page's regions, found again on the volume read
    #      back: the clamp and the lower half of the bore.
    encastrement = pc.mesh.points_on_plane(volume, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0])
    bas_axe = [LENGTH, -THICKNESS, HEIGHT / 2.0]
    haut_axe = [LENGTH, 2.0 * THICKNESS, HEIGHT / 2.0]
    noeuds_alesage = pc.mesh.points_on_cylinder(peau, bas_axe, haut_axe, HOLE_RADIUS)
    alesage = pc.mesh.elements_on(peau, noeuds_alesage, strict=True)
    z_alesage = pc.node_field.positions(alesage, ["Z"])
    noeuds_bas = pc.mesh.select(z_alesage, le=HEIGHT / 2.0 + TOL)
    demi_alesage = pc.mesh.consolidate(
        pc.mesh.elements_on(alesage, noeuds_bas, strict=True)
    )

    # FR — Le modèle de la page « Élasticité » ; le matériau porte en plus
    #      « alpha », que la dilatation lira.
    # EN — The "Élasticité" page's model; the material also carries "alpha",
    #      which the expansion will read.
    fes = pc.FiniteElementSpace(volume)
    elasticite = pc.model.elasticity(fes, "full_3d")
    blocage_x = pc.model.dirichlet(encastrement, elasticite, "u_x")
    blocage_y = pc.model.dirichlet(encastrement, elasticite, "u_y")
    blocage_z = pc.model.dirichlet(encastrement, elasticite, "u_z")
    modele = elasticite | blocage_x | blocage_y | blocage_z
    materiaux = pc.element_field.material_field(
        modele, [("E", E), ("nu", NU), ("alpha", ALPHA)]
    )
    K = pc.matrix.stiffness(modele, materiaux)

    demi_fes = pc.FiniteElementSpace(demi_alesage)
    verticale = pc.model.flux(demi_fes, elasticite, "f_z")
    densite_verticale = pc.element_field.material_field(
        verticale, [("phi_f_z", -PRESSION)]
    )
    force_surfacique = pc.node_field.external_forces(verticale, densite_verticale)

    # FR — La température aux points de Gauss, puis la dilatation
    #      ε_th = α (T − T_ref), la contrainte qu'elle donnerait si rien ne la
    #      gênait, [C] ε_th, et ses forces nodales ∫ [B]ᵀ [C] ε_th dV.
    # EN — The temperature at the Gauss points, then the expansion
    #      ε_th = α (T − T_ref), the stress it would give were nothing in its
    #      way, [C] ε_th, and its nodal forces ∫ [B]ᵀ [C] ε_th dV.
    t_gauss = pc.element_field.interp_to_gauss(temperature, fes)
    eps_th = pc.element_field.thermal_strain(t_gauss, materiaux, fes, T_REF)
    sigma_th = pc.element_field.elastic_stress(modele, eps_th, materiaux)
    force_th = (
        pc.node_field.divergence(sigma_th, "sigma")
        .rename_component("div_sigma_x", "f_x")
        .rename_component("div_sigma_y", "f_y")
        .rename_component("div_sigma_z", "f_z")
    )

    # FR — Les deux chargements partagent les nœuds de l'alésage : la force
    #      surfacique est ramenée sur le support de la dilatation, puis `+`
    #      somme nœud à nœud.
    # EN — Both loads share the bore's nodes: the surface force is brought
    #      onto the expansion's support, then `+` sums node by node.
    second_membre = force_th + pc.node_field.restrict_like(force_surfacique, force_th)

    u = pc.solver.solve(K, second_membre)
    r_x = -pc.node_field.restrict(u, blocage_x.multiplier_mesh()).sum("lambda")
    r_y = -pc.node_field.restrict(u, blocage_y.multiplier_mesh()).sum("lambda")
    r_z = -pc.node_field.restrict(u, blocage_z.multiplier_mesh()).sum("lambda")

    coords = volume.coords()
    deformee = coords.add_config("déformée ×50")
    pc.coords.displace_from(u * AMPLIFICATION, 0, deformee)
    with coords.in_config(deformee):
        volume.plot(
            view=VUE,
            title="Dilatation et force surfacique — déformée ×50, u_x (mm)",
            field=u * 1000.0,
            component="u_x",
            cmap="viridis",
            save=os.path.join(OUT, "mecanique-dilatation-deformee.svg")
            if OUT
            else None,
        )

    # FR — Seule la part de la déformation que la dilatation n'explique pas
    #      crée une contrainte : σ = [C] (ε − ε_th).
    # EN — Only the part of the strain the expansion does not account for
    #      creates a stress: σ = [C] (ε − ε_th).
    deformation = pc.element_field.strain(modele, u, materiaux)
    contrainte = pc.element_field.elastic_stress(
        modele, deformation - eps_th, materiaux
    )
    von_mises_mpa = pc.element_field.von_mises(modele, contrainte) * 1e-6
    volume.plot(
        view=VUE,
        title="Dilatation et force surfacique — von Mises (MPa)",
        field=von_mises_mpa,
        component="von_mises",
        cmap="viridis",
        save=os.path.join(OUT, "mecanique-dilatation-von-mises.svg") if OUT else None,
    )

    print(f"relu : {volume.cell_count} HEX20, {len(evolution.frames())} champs")
    print(
        f"température relue : de {temperature.min('T'):.2f} à "
        f"{temperature.max('T'):.2f} °C"
    )
    print(
        f"forces de la dilatation, sommées : {force_th.sum('f_x'):.2e} "
        f"{force_th.sum('f_y'):.2e} {force_th.sum('f_z'):.2e} N"
    )
    print(f"réactions : {r_x:.1f} {r_y:.1f} {r_z:.1f} N")
    print(
        f"u_x de {u.min('u_x') * 1000.0:.3f} à {u.max('u_x') * 1000.0:.3f} mm, "
        f"u_z de {u.min('u_z') * 1000.0:.3f} à {u.max('u_z') * 1000.0:.3f} mm"
    )
    print(
        f"aux points de Gauss : von Mises de {von_mises_mpa.min('von_mises'):.2f} "
        f"à {von_mises_mpa.max('von_mises'):.2f} MPa"
    )


if __name__ == "__main__":
    main()

Suite : Dilatation, matériau hétérogène.