3 Application numérique
3.1 Importation des librairies
Cette section organise les imports par catégorie pour plus de lisibilité :
- Bibliothèques standard :
math,timepour les opérations de base - Calcul scientifique :
numpypour le calcul numérique et les tableaux - FEniCSx : Framework Elements Finis pour la résolution du problème mécanique
- UFL (Unified Form Language) : Langage pour définir les formes variationnelles
- Gmsh : Générateur de maillage pour créer la géométrie du cylindre
- Visualisation :
matplotlib,pyvista,panelpour les graphiques et rendu 3D - Calcul parallèle :
mpi4py,petsc4pypour la parallélisation FEniCS
3.2 Géométrie et maillage du cylindre
3.2.1 Création de la géométrie avec Gmsh
La géométrie est définie par un cylindre elliptique de demies axes \(\rho_x\) et \(\rho_y\) et de hauteur \(h\) :
- On crée d’abord une ellipse dans le plan \(xy\) avec l’équation \(\frac{x^2}{\rho_x^2} + \frac{y^2}{\rho_y^2} = 1\)
- Puis on extrude cette surface selon l’axe \(z\) pour obtenir un cylindre elliptique 3D
3.2.2 Maillage
- Le paramètre
lccontrôle la finesse du maillage (Characteristic Length) - Un maillage plus fin donne des résultats plus précis mais augmente le temps de calcul
- Le maillage généré est ensuite converti au format XDMF pour FEniCSx
3.3 Définition de l’espace fonctionnel
3.3.1 Espaces d’approximation
En éléments finis, l’espace fonctionnel \(V\) est un espace de fonctions sur lequel on cherche la solution. Pour un problème 3D en déplacement :
\[V = \{ \mathbf{u} \in (H^1(\Omega))^3 \}\]
- \(CG\), 1 : Éléments de type Lagrange continu d’ordre 1 (linéaires par tétraèdres)
- \((domain.geometry.dim,)\) : Vecteur 3D (une composante par direction \(x, y, z\))
3.3.2 Fonctions d’essai et de test
- TrialFunction \(u\) : Représente le déplacement inconnu \(\mathbf{u}\)
- TestFunction \(v\) : Fonction test pour la formulation variationnelle
Ces fonctions appartiennent à l’espace \(V\) et satisfont les conditions aux limites.
V = fem.functionspace(domain, ("CG", 1, (domain.geometry.dim,)))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)3.4 Définition des frontières du domaine
3.4.1 Identification géométrique des surfaces
Pour appliquer les conditions aux limites, on identifie les surfaces du cylindre :
- Face inférieure \(\Sigma_0\) : \(z = 0\) (base du cylindre)
- Face supérieure \(\Sigma_h\) : \(z = h\) (sommet du cylindre)
- Surface latérale \(\Sigma_l\) : Surface elliptique \(\frac{x^2}{\rho_x^2} + \frac{y^2}{\rho_y^2} = 1\)
3.5 Visualisation des frontières du domaine
3.5.1 Validation du maillage
Cette visualisation permet de vérifier que :
- Les marqueurs de frontières sont correctement attribués
- Chaque surface est bien identifiée (inférieure, supérieure, latérale)
- Le maillage est régulier et sans défaut
Les couleurs permettent de distinguer les différentes surfaces géométriques qui seront utilisées pour les conditions aux limites.
3.6 Conditions aux limites
3.6.1 Conditions de Dirichlet (encastrement)
Une condition aux limites de Dirichlet impose la valeur de la solution sur une partie du bord :
\[\mathbf{u} = \mathbf{u}_D \quad \text{sur } \Gamma_D\]
Dans ce problème, on applique un encastrement (déplacement nul) sur les surfaces sélectionnées :
- \(\mathbf{u}_D = (0, 0, 0)\) : déplacement nul dans toutes les directions
3.6.2 Localisation des degrés de liberté
fem.locate_dofs_topological(V, fdim, facets) identifie les degrés de liberté (DOFs) correspondant aux facettes géométriques marquées.
# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
# Application des conditions aux limites sur chaque face
bc0_S0 = fem.dirichletbc(
u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V
) # Face inférieure
bc0_Sh = fem.dirichletbc(
u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V
) # Face supérieure
bc0_Sl = fem.dirichletbc(
u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V
) # Surface latérale3.7 Propriétés matériau
3.7.1 Coefficients de Lamé
Pour un matériau élastique linéaire isotrope, le comportement est caractérisé par les coefficients de Lamé :
- \(\mu\) (module de cisaillement) : \(\mu = \frac{E}{2(1+\nu)}\)
- \(\lambda\) (premier paramètre de Lamé) : \(\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}\)
où \(E\) = Module de Young (ici \(E = 210\) GPa pour l’acier) et \(\nu\) = Coefficient de Poisson (ici \(\nu = 0.3\))
3.7.2 Loi de Hooke généralisée (3D)
La relation contrainte-déformation en 3D s’écrit :
\[\boldsymbol{\sigma} = \lambda \, \text{tr}(\boldsymbol{\varepsilon}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}\]
Code
def calculate_lame_constants(E, nu):
"""
Calcule les coefficients de Lamé (mu et lambda) à partir de E et nu.
Args:
E (float): Module de Young (Pa).
nu (float): Coefficient de Poisson (sans unité).
Returns:
tuple: (mu, lambda) en Pascals.
"""
mu = E / (2 * (1 + nu))
lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
return mu, lambda_mu, lambda_ = display_elastic_constants("value")
# Formater pour afficher en GPa
print(f"Module de cisaillement (μ) = {mu / 1e9:.2f} GPa")
print(f"Premier paramètre de Lamé (λ) = {lambda_ / 1e9:.2f} GPa")Module de cisaillement (μ) = 80.77 GPa
Premier paramètre de Lamé (λ) = 121.15 GPa3.8 Définition des chargements
3.8.1 Forces volumiques (pesanteur)
Une force volumique \(\mathbf{f}\) agit sur tout le volume du matériau. Pour la pesanteur :
\[\mathbf{f} = \rho \, \mathbf{g}\]
où \(\rho\) = masse volumique (ici \(\rho = 7850\) kg/m³ pour l’acier) et \(\mathbf{g}\) = accélération de la gravité (\(g = 9.81\) m/s²)
On applique une force gravitationnelle selon l’axe \(z\) négatif : \(\mathbf{f} = (0, 0, -\rho g)\)
3.8.2 Tenseur des déformations
Le tenseur de déformation linéarisé est défini par :
\[\boldsymbol{\varepsilon}(\mathbf{u}) = \text{sym}(\nabla\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)\]
Code
def calculate_force_volumique(rho, g):
"""
Calcule la force volumique f à partir de la masse volumique rho
et de l'accélération gravitationnelle g.
"""
return rho * g
# données
rho = 7850 # Masse volumique de l'acier en kg/m^3
g = 9.81 # Accélération gravitationnelle en m/s^2
fv = calculate_force_volumique(rho, g)
# Définition de la force volumique (nulle dans ce cas)
f = fem.Constant(domain, default_scalar_type((0, 0, -fv)))3.9 Fonctions d’essai et de test
3.9.1 Approximation variationnelle
Dans la méthode des éléments finis, on cherche une approximation \(\mathbf{u}_h\) de la solution exacte dans un espace de dimension finie :
\[\mathbf{u}_h = \sum_{i=1}^{N} u_i \boldsymbol{\phi}_i\]
où \(\boldsymbol{\phi}_i\) sont les fonctions de forme et \(u_i\) les inconnues (coefficients à déterminer).
- TrialFunction \(u\) : Représente la solution inconnue (les coefficients \(u_i\))
- TestFunction \(v\) : Fonction test qui multiplie l’équation différentielle pour obtenir la formulation variationnelle
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)3.10 Forme bilinéaire \(a\)
3.10.1 Énergie de déformation
La forme bilinéaire \(a(\mathbf{u}, \mathbf{v})\) représente l’opérateur rigide (matrice de rigidité) du problème élastique :
\[a(\mathbf{u}, \mathbf{v}) = \int_\Omega \boldsymbol{\sigma}(\mathbf{u}) : \boldsymbol{\varepsilon}(\mathbf{v}) \, d\Omega\]
En utilisant la loi de Hooke, cela devient :
\[a(\mathbf{u}, \mathbf{v}) = \int_\Omega \left[ \lambda \, \text{tr}(\boldsymbol{\varepsilon}(\mathbf{u})) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}) \right] : \boldsymbol{\varepsilon}(\mathbf{v}) \, d\Omega\]
Cette forme est bilinéaire (linéaire en \(\mathbf{u}\) et en \(\mathbf{v}\)) et symétrique.
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx # Forme bilinéaire3.11 Forme linéaire \(L\)
3.11.1 Travail des forces extérieures
La forme linéaire \(L(\mathbf{v})\) représente le second membre (terme source) du problème variationnel :
\[L(\mathbf{v}) = \int_\Omega \mathbf{f} \cdot \mathbf{v} \, d\Omega + \int_{\Gamma_N} \mathbf{T} \cdot \mathbf{v} \, d\Gamma\]
où \(\mathbf{f}\) = force volumique (pesanteur) et \(\mathbf{T}\) = force surfacique (traction sur \(\Gamma_N\))
Dans ce problème, on considère uniquement la force volumique due à la gravité :
\[L(\mathbf{v}) = \int_\Omega (0, 0, -\rho g) \cdot \mathbf{v} \, d\Omega\]
Cette forme est linéaire en \(\mathbf{v}\).
Code
def update_linear_form(change):
global L
selected_surface = change["new"]
# Mise à jour de la forme linéaire
L = ufl.dot(f, v) * ufl.dx
# Mise à jour de l'affichage
with output:
output.clear_output(wait=True)
print(f"force gravitationnelle exercée sur {selected_surface}")3.12 Résolution par la formule variationnelle
3.12.1 Problème variationnel
Le problème aux éléments finis s’écrit :
Trouver \(\mathbf{u} \in V\) tel que :
\[a(\mathbf{u}, \mathbf{v}) = L(\mathbf{v}) \quad \forall \mathbf{v} \in V_0\]
où \(V_0\) est l’espace des fonctions tests satisfaisant les conditions aux limites de Dirichlet homogènes.
3.12.2 Méthode de résolution
- LinearProblem : Assemble le système linéaire \(KU = F\) à partir des formes \(a\) et \(L\)
- KSP CG : Solveur itératif (Conjugate Gradient) adapté aux matrices symétriques définies positives
- Préconditionneur Jacobi : Préconditionneur diagonal pour accélérer la convergence
Code
# Fonction pour mettre à jour le problème
def update_problem(change):
global problem
selected_bc = change["new"]
problem = LinearProblem(
a,
L,
bcs=bc_dict[selected_bc],
petsc_options={"ksp_type": "cg", "pc_type": "jacobi"},
)
# Mise à jour de l'affichage
with output:
output.clear_output(wait=True)
print(f"Encastrement appliqué sur {selected_bc}")3.13 Résolution du problème
uh = problem.solve()Temps de calcul pour la résolution : 0.1908 secondes3.14 Visualisation des déplacements
3.14.1 Interprétation des résultats
La solution numérique \(\mathbf{u}_h\) donne le champ de déplacement en chaque nœud du maillage :
\[\mathbf{u}_h = (u_x, u_y, u_z)\]
Pour visualiser les résultats :
- Warping : On déforme géométriquement le maillage par le champ de déplacement (facteur d’exaggération)
- Coloration : Les couleurs représentent l’amplitude du déplacement (ici \(u_z\))
3.14.2 Analyse physique
- Sous l’effet du poids propre (gravité), le cylindre se comprime axialement
- Le déplacement \(u_z\) est négatif (vers le bas) et maximal au sommet libre
- La base encastrée (\(u_z = 0\)) ne bouge pas
- Forme parabolique : La déformée \(u_z(z)\) suit une loi parabolique car l’effort axial varie linéairement avec \(z\) (nul au sommet libre, maximal à la base encastrée), ce qui donne \(u_z(z) \propto z^2\)
3.15 Visualisation des contraintes de von Mises
3.15.1 Critère de von Mises
Le critère de von Mises prédit la limite élastique des matériaux ductiles :
\[\sigma_{vm} = \sqrt{\frac{1}{2} \left[ (\sigma_{xx}-\sigma_{yy})^2 + (\sigma_{yy}-\sigma_{zz})^2 + (\sigma_{zz}-\sigma_{xx})^2 + 6(\sigma_{xy}^2 + \sigma_{yz}^2 + \sigma_{zx}^2) \right]}\]
En notation tensorielle :
\[\sigma_{vm} = \sqrt{3J_2}\]
où \(J_2\) est le deuxième invariant du tenseur déviateur des contraintes.
3.15.2 Interprétation
- \(\sigma_{vm} < \sigma_e\) : Le matériau reste élastique
- \(\sigma_{vm} \geq \sigma_e\) : Risque de plastification
3.16 Comparaison Analytique-Numérique
3.16.1 Validation du modèle numérique
Pour valider les résultats FEM, on compare avec la solution analytique du déplacement axial d’un cylindre sous son poids propre :
\[u_z^{théorique} = -\frac{\rho g h^2}{2E}\]
3.16.2 Analyse de convergence
L’erreur relative est calculée par :
\[e_{rel} = \frac{|u_z^{théorique} - u_z^{EF}|}{|u_z^{théorique}|} \times 100\%\]
On observe que : - L’erreur diminue progressivement avec le raffinement du maillage - La convergence est typique d’une méthode d’éléments finis d’ordre 1 - Un maillage plus fin donne des résultats plus proches de la solution analytique
