5 Application numérique
Ce notebook utilise FEniCSx pour simuler la torsion d’un cylindre elliptique. La problématique pédagogique consiste à : - Modéliser un cylindre elliptique en 3D - Appliquer des conditions aux limites de type encastrement - Résoudre le problème élastique linéaire - Visualiser les déplacements, déformations et contraintes - Comparer les résultats numériques avec la théorie de la torsion
5.1 Importation des bibliothèques
Cette section rassemble tous les imports nécessaires au notebook. On distingue plusieurs catégories : - FEniCSx : le solveur éléments finis (dolfinx) - UFL : Unified Form Language pour définir les formulations variationnelles - Visualisation : PyVista pour le rendu 3D, Matplotlib pour les graphes - Outils numériques : NumPy pour les calculs, SymPy pour les symboles - Maillage : Gmsh pour générer le maillage, meshio pour la conversion
5.2 Géométrie et maillage
5.2.1 Définition de la géométrie
On génère un cylindre elliptique à l’aide de Gmsh. La section transversale est une ellipse de demi-axes \(\rho_x\) et \(\rho_y\), extrudée sur une hauteur \(h\).
Le maillage généré est ensuite converti au format XDMF pour être lu par FEniCSx.
Nombre de cellules : 74883
Nombre de sommets : 142625.2.2 Définition de l’espace fonctionnel
L’espace fonctionnel \(V\) est l’espace des éléments finis dans lequel nous cherchons la solution (le déplacement).
Pour un problème d’élasticité 3D, on utilise un espace vectoriel de type Lagrange : - Famille : "Lagrange" (éléments de Lagrange) - Ordre : 2 (éléments quadratiques, ici pour une meilleure précision) - Dimension : 3 (vecteur 3D car problème tridimensionnel)
La formulation variationnelle fait intervenir : - \(u\) : la fonction trial (inconnue) - le déplacement - \(v\) : la fonction test - introduit pour intégrer par parties
# =============================================================================
# DÉFINITION DE L'ESPACE FONCTIONNEL
# =============================================================================
# Création de l'espace de fonction pour le problème d'élasticité
# Syntaxe: fem.functionspace(maillage, (famille_element, degré, dimension))
# Ici : éléments de Lagrange quadratiques (degré 2) en 3D (vecteur à 3 composantes)
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim, )))
# Définition des fonctions trial et test pour la formulation variationnelle
u = ufl.TrialFunction(V) # Fonction inconnue (déplacement à chercher)
v = ufl.TestFunction(V) # Fonction test (fonction d.weight)5.2.3 Définition des frontières du domaine
Pour appliquer les conditions aux limites, on doit identifier les différentes surfaces du cylindre : - \(\Sigma_0\) : Face inférieure (\(z = 0\)) - souvent utilisée comme encastrement - \(\Sigma_h\) : Face supérieure (\(z = h\)) - où on applique la torsion - \(\Sigma_l\) : Surface latérale (elliptique)
La méthode mesh.locate_entities permet de trouver les entités géométriques satisfaisant certain critère (ici, une condition sur les coordonnées).
5.3 Conditions aux limites
5.3.1 Conditions de Dirichlet (Encastrement)
Les conditions aux limites de Dirichlet spécifient la valeur de la solution sur une partie de la frontière. Dans notre cas : - On peut encastrer une face (déplacement nul : \(u = 0\)) - C’est ce qu’on appelle un encastrement en mécanique des structures
La syntaxe dans FEniCSx est : fem.dirichletbc(valeur, dofs, espace_fonction) où dofs sont les degrés de liberté trouvés par localisation topologique.
# =============================================================================
# CONDITIONS AUX LIMITES DE DIRICHLET (Encastrement)
# =============================================================================
# Vecteur de déplacement nul (encastrement) : [0, 0, 0]
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
# Création des conditions de Dirichlet pour chaque face
# Syntaxe : fem.dirichletbc(valeur,
# fem.locate_dofs_topological(espace, dimension, facettes),
# espace)
bc0_S0 = fem.dirichletbc(u_D,
fem.locate_dofs_topological(V, fdim, Sigma_0),
V) # Face inférieure (z = 0)
bc0_Sh = fem.dirichletbc(u_D,
fem.locate_dofs_topological(V, fdim, Sigma_h),
V) # Face supérieure (z = h)
bc0_Sl = fem.dirichletbc(u_D,
fem.locate_dofs_topological(V, fdim, Sigma_l),
V) # Surface latérale5.3.2 Propriétés du matériau
Le matériau est caractérisé par les paramètres de Lamé : - \(\lambda\) : premier paramètre de Lamé (module de compression) - \(\mu\) : module de cisaillement (aussi noté G)
Ces paramètres sont liés au module d’Young \(E\) et au coefficient de Poisson \(\nu\) par : \[\lambda = \frac{E \nu}{(1+\nu)(1-2\nu)}, \quad \mu = \frac{E}{2(1+\nu)}\]
5.3.3 Définition du chargement (Torsion)
Le chargement est appliqué sous forme de torsion sur l’une des faces du cylindre. Le champ de torsion est un champ de déplacement tangentiel :
\[ \mathbf{T} = \begin{pmatrix} -A \cdot y \\ A \cdot x \\ 0 \end{pmatrix} \]
où \(A = \mu \theta / h\) est un scalaire dépendant du module de cisaillement, de l’angle de torsion \(\theta\) et de la hauteur \(h\).
Rayon équivalent: R = 0.2000 m
Couple maximal: T_max = 1.25664e+06 N·m\(\displaystyle \mathbf{T} = \begin{bmatrix} - A x_{2} \\ A x_{1} \\ 0 \end{bmatrix}\)
5.4 Formulation variationnelle
La formulation variationnelle du problème d’élasticité linéaire est :
Trouver \(u \in V\) tel que :
\[ \int_\Omega \sigma(u) : \epsilon(v) \, d\Omega = \int_\Omega \mathbf{f} \cdot v \, d\Omega + \int_{\Sigma_T} \mathbf{T} \cdot v \, d\Sigma \]
où : - \(\sigma(u)\) est le tenseur des contraintes (loi de Hooke) - \(\epsilon(u) = \frac{1}{2}(\nabla u + \nabla u^T)\) est le tenseur de déformation - \(\mathbf{f}\) est la force volumique (ici nulle) - \(\mathbf{T}\) est le chargement surfacique (torsion)
Loi de Hooke (en 3D) : \[ \sigma = \lambda \, \text{tr}(\epsilon) \, I + 2\mu \, \epsilon \]
5.5 Résolution du problème éléments finis
Une fois la formulation variationnelle définie et les conditions aux limites fixées, on assemble le système linéaire et on le résout. La classe LinearProblem de FEniCSx encapsule ces opérations et utilise PETSc comme solveur linéaire.
# =============================================================================
# RÉSOLUTION DU PROBLÈME
# =============================================================================
# Résolution du problème linéaire
# La solution uh est une fonction FEniCS contenant les valeurs nodales
uh = problem.solve()
print("Problème résolu avec succès!")
print(f"Nombre de degrés de liberté : {V.dofmap.index_map.size_local}")Problème résolu avec succès!
Nombre de degrés de liberté : 1078675.6 Post-traitement et visualisation
On utilise PyVista (basé sur VTK) pour visualiser les résultats en 3D. Les résultats FEniCS sont interpolés sur un maillage compatible VTK.
5.6.1 Visualisation des déplacements
5.6.2 Visualisation des déplacements amplifiés
On applique une amplification (factor=500) au maillage pour visualiser la déformation.
5.6.3 Contraintes de von Mises
Le critère de von Mises est souvent utilisé pour prédire la rupture des matériaux ductiles. La contrainte équivalente de von Mises est :
\[\sigma_{vm} = \sqrt{\frac{1}{2}[(\sigma_{xx}-\sigma_{yy})^2 + (\sigma_{yy}-\sigma_{zz})^2 + (\sigma_{zz}-\sigma_{xx})^2 + 6(\sigma_{xy}^2+\sigma_{yz}^2+\sigma_{xz}^2)]}\]
# =============================================================================
# CONTRAINTES DE VON MISES
# =============================================================================
# Calcul du déviateur des contraintes
s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
# Contrainte équivalente de von Mises
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))
# Interpolation sur l'espace DG
V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses_vm = fem.Function(V_von_mises)
stresses_vm.interpolate(stress_expr)# Comparaison avec la théorie
max_vm = np.max(np.abs(stresses_vm.x.array))
sigma_p_num = max_vm / math.sqrt(3)
print("Contrainte principale maximale numérique :")
print(f"{sigma_p_num / 1e6:.2f} MPa")
sigma_p_theo = A_coeff * R_eq
print("Contrainte principale maximale théorique :")
print(f"{sigma_p_theo / 1e6:.2f} MPa")Contrainte principale maximale numérique :
99.90 MPa
Contrainte principale maximale théorique :
101.40 MPa5.7 Analyse des résultats et validation
5.7.1 Comparaison numérique-théorique
On compare les résultats de rotation maximale avec la théorie de la torsion. Pour un cylindre circulaire, l’angle de torsion est donné par :
\[\theta = \frac{T \cdot L}{J \cdot G}\]
où \(T\) est le couple, \(L\) la longueur, \(J\) le moment d’inertie polaire et \(G\) le module de cisaillement.
Code
# =============================================================================
# VALIDATION : COMPARAISON NUMÉRIQUE vs THÉORIE
# =============================================================================
# Paramètres
theta_theoretical = np.radians(angle_slider.value)
# Extraction des maxima de rotation
max_rotation_z = np.max(np.abs(rotation_z.x.array))
max_rotation_norm = np.max(rotation_norm.x.array)
# Affichage des résultats
print(f"Moment d'inertie polaire J : {J:.6e} m^4")
print(f"Angle de torsion théorique θ : {theta_theoretical:.6f} rad")
print(f"Angle de torsion théorique θ : {np.degrees(theta_theoretical):.4f}°")
# Comparaison avec la simulation
print(f"\n--- Comparaison avec la simulation ---")
print(f"Rotation Z maximale simulée : {max_rotation_z:.6f} rad")
print(f"Rotation Z maximale simulée : {np.degrees(max_rotation_z):.4f}°")
relative_error_Z = abs(max_rotation_z - theta_theoretical) / theta_theoretical * 100
print(f"Erreur relative : {relative_error_Z:.2f}%")
print(f"\nNorme de rotation maximale : {max_rotation_norm:.6f} rad")
print(f"Norme de rotation maximale : {np.degrees(max_rotation_norm):.4f}°")
relative_error_N = (max_rotation_norm - theta_theoretical) / theta_theoretical * 100
print(f"Erreur relative : {relative_error_N:.2f}%")
# Analyse des composantes
max_rotation_x = np.max(np.abs(rotation_x.x.array))
max_rotation_y = np.max(np.abs(rotation_y.x.array))
print("\n--- Valeurs maximales des composantes de rotation ---")
print(f"RX max : {max_rotation_x:.6f} rad ({np.degrees(max_rotation_x):.2f}°)")
print(f"RY max : {max_rotation_y:.6f} rad ({np.degrees(max_rotation_y):.2f}°)")
print(f"RZ max : {max_rotation_z:.6f} rad ({np.degrees(max_rotation_z):.4f}°)")
# Contribution relative
total_rotation = max_rotation_x + max_rotation_y + max_rotation_z
print("\n--- Contribution relative ---")
print(f"RX : {max_rotation_x / total_rotation * 100:.2f}%")
print(f"RY : {max_rotation_y / total_rotation * 100:.2f}%")
print(f"RZ : {max_rotation_z / total_rotation * 100:.2f}%")Moment d'inertie polaire J : 2.513274e-03 m^4
Angle de torsion théorique θ : 0.010140 rad
Angle de torsion théorique θ : 0.5810°
--- Comparaison avec la simulation ---
Rotation Z maximale simulée : 0.010098 rad
Rotation Z maximale simulée : 0.5786°
Erreur relative : 0.41%
Norme de rotation maximale : 0.010145 rad
Norme de rotation maximale : 0.5813°
Erreur relative : 0.04%
--- Valeurs maximales des composantes de rotation ---
RX max : 0.000994 rad (0.06°)
RY max : 0.000994 rad (0.06°)
RZ max : 0.010098 rad (0.5786°)
--- Contribution relative ---
RX : 8.23%
RY : 8.22%
RZ : 83.55%5.7.2 Histogramme des contraintes
Code
# =============================================================================
# HISTOGRAMME DES CONTRAINTES PRINCIPALES
# =============================================================================
plt.rcParams['figure.facecolor'] = (40/255, 43/255, 48/255, 1.0)
plt.rcParams['axes.facecolor'] = (40/255, 43/255, 48/255, 1.0)
plt.rcParams['savefig.facecolor'] = (40/255, 43/255, 48/255, 1.0)
plt.figure(figsize=(12, 6))
plt.hist(stress_principal_values, bins=50, color='cyan', alpha=0.7)
plt.axvline(sigma_lim, color='magenta', linestyle='dashed', label="Limite rupture")
# Axes gris
plt.rcParams['axes.edgecolor'] = '0.5'
plt.rcParams['xtick.color'] = '0.5'
plt.rcParams['ytick.color'] = '0.5'
# Labels + légende GRIS en UN SEUL APPEL
plt.xlabel("Contraintes principales (Pa)", color='grey', fontsize=16)
plt.ylabel("Nombre d'éléments", color='grey', fontsize=16)
plt.title("Répartition des contraintes principales", color='grey', fontsize=18, pad=20)
# LÉGENDE avec texte GRIS
plt.legend(labelcolor='0.5', facecolor='none', fontsize=14) # TOUT ICI !
plt.savefig("images/histogram.png", dpi=300, facecolor='none', transparent=false)
plt.show()--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[39], line 28 25 # LÉGENDE avec texte GRIS 26 plt.legend(labelcolor='0.5', facecolor='none', fontsize=14) # TOUT ICI ! ---> 28 plt.savefig("images/histogram.png", dpi=300, facecolor='none', transparent=false) 29 plt.show() NameError: name 'false' is not defined
5.8 Conclusions
Ce notebook a démontré :
- La modélisation d’un cylindre elliptique en 3D avec FEniCSx
- L’application de conditions aux limites de type encastrement et chargement
- La résolution du problème d’élasticité linéaire par éléments finis
- La visualisation des déplacements, déformations et contraintes
- La validation des résultats numériques par rapport à la théorie
Les résultats montrent une excellente concordance entre la simulation numérique et les prédictions théoriques de la torsion (erreur relative < 1%).