Éléments théoriques de base

Généralités

L'optimisation topologique consiste à trouver la répartition optimale de la matière dans une pièce donnée soumise à des chargements et sous certaines limitations. L'inconnue de ce type de problème est la topologie, c'est-à-dire comment la matière est-elle répartie ?

Dans un cadre mécanique, ce type de problème peut être formulé en : trouver la distribution optimale de rigidité \(\mathbfcal{k}\), ce qui peut s'écrire :

(1)\[\begin{split}\min_{\mathbfcal{k}(m)} \quad & f(\mathbfcal{k}) & \\ \textsf{tel que} \quad & g_i(\mathbfcal{k}) &= 0 \\\end{split}\]

Explicitons certains termes :

  • Par répartition optimale de la matière on entend la distribution spatiale de rigidité \(\mathbfcal{k}(m)\) en tout point \(m \in \Omega\) qui minimise la fonction objectif \(f\), indicatrice du souhait de l'utilisateur.

  • Par pièce donnée on entend que la recherche de la topologie se fait sur un domaine limité de l'espace \(\Omega\).

  • Par chargements on entend que cette pièce subit des conditions aux limites, par exemple des déplacements imposés et/ou des forces appliquées sur certaines zones.

  • Par limitations on entend que l'optimisation se fait sous contraintes, représentées par l'ensemble de fonctions \(g_i(\mathbfcal{k}) = 0\).

Le problème est illustré ci-dessous dans le cas d'une poutre en flexion 3 points.

../_images/99_lines_img_1.png

Problème d'optimisation topologique sur une poutre en flexion [SIGMUND-2001]

Plusieurs choix de fonctions objectif \(f\) et contraintes \(g_i\) sont alors possibles pour optimiser une structure, comme par exemple :

  • chercher la pièce la plus rigide possible, sans dépasser un certain volume ;

  • chercher la pièce la plus légère possible, sans dépasser une certaine fréquence de résonnance ;

  • chercher la pièce avec les contraintes mécaniques les plus faibles, sans dépasser un certain niveau de déplacement.

Un choix très répendu est de minimiser la compliance (c'est-à-dire l'énergie de déformation ou encore le travail des forces extérieures) sous contrainte de volume.

Dans ce document nous illustrerons brièvement la méthode SIMP (pour Solid Isotropic Material with Penalization) qui est la méthode d'optimisation topologique la plus répendue dans les codes industriels et celle mise en oeuvre dans Cast3M via la procédure TOPOPTIM. Le lecteur intéressé pourra consulter de nombreux ouvrages sur le sujet, comme par exemple :

  • Les livres de référence [BENDSOE-1995] et [BENDSOE-AL-2004] qui détaillent rigouresement la théorie derrière l'optimisation topologique.

  • Les articles pédagogiques [SIGMUND-2001] et [ANDREASSEN-AL-2011] qui présentent des implémentations sur Matlab, en respectivement 99 lignes puis 88 lignes, d'un algorithme d'optimisation topologique. La procédure TOPOPTIM de Cast3M en est inspiré et l'exemple utilisé dans ce document en est une reproduction.

La méthode SIMP

Le domaine de conception est discrétisé selon un maillage de \(N\) éléments. Chaque élément \(e\) se voit attribué une variable de densité \(x_e\) et un module d'Young \(E_e\) selon une loi puissance :

(2)\[E_e(x_e) = E_{\textrm{min}} + {x_e}^p(E_0-E_{\textrm{min}})\]

avec :

  • \(x_e \in [0,1]\) la densité de l'élément \(e\) (où \(x_e=0\) correspond à un élément "vide" et \(x_e=1\) un élément "plein")

  • \(E_0\) le module d'Young du matériau

  • \(E_{\textrm{min}}\) une valeur "faible" de module d'Young choisie pour éviter que la matrice de rigidité du modèle éléments finis devienne singulière (typiquement \(E_{\textrm{min}}=E_0.10^{-6}\))

  • \(p\) un coefficient de pénalisation (typiquement \(p=3\)) introduit afin de défavoriser les valeurs intermédiaires \(0<x_e<1\) lors de l'optimisation

Le problème d'optimisation topologique de minimisation de la compliance, sous contrainte de volume, s'écrit alors :

(3)\[\begin{split}\min_{\textbf{x}} \quad & f(\textbf{x}) = \textbf{U}^T.\mathbfcal{K}.\textbf{U} = \sum_{e=1}^{N} E_e(x_e)\textbf{u}_e^T.\mathbfcal{k}_0.\textbf{u}_e \\ \textsf{tel que} \quad & g(\textbf{x}) = \phi(\textbf{x}) - \phi_0 = \sum_{e=1}^N \frac{x_eV_e}{V_0} - \phi_0 = 0 \\ & \mathbfcal{K}.\textbf{U} = \textbf{F} \\ & 0 \leq x_e \leq 1 \quad \forall e \in \{1;\dots;N\} \\\end{split}\]

Ce problème de minimisation de la compliance est équivalent à maximiser la raideur de la structure. Précisons quelques variables :

  • \(\textbf{x}\) est le vecteur des densités des \(N\) éléments du maillage, les iconnues du problème, encore appellées variables de conception

  • \(f\) est la fonction objectif, la compliance de la structure (homogène à une énergie)

  • \(\textbf{U}\) et \(\textbf{F}\) sont les vecteurs déplacements et forces aux noeuds du maillage

  • \(\textbf{u}_e\) est le vecteur déplacements réduit à l'élément \(e\)

  • \(\mathbfcal{K}=\sum_{e=1}^N E(x_e)\mathbfcal{k}_0\) est la matrice de rigidité globale, assemblée sur les \(N\) éléments du maillage

  • \(\mathbfcal{k}_0\) est la matrice de rigidité élémentaire pour un module d'Young unitaire

  • \(g(\textbf{x})\) est la fonction contrainte d'égalité

  • \(\phi(\textbf{x})\) est la fraction volumique de la topologie et \(\phi_0\) est la fraction volumique imposée

  • \(V_e\) est le volume de l'élément \(e\) et \(V_0\) est le volume du domaine de conception \(\Omega\)

Résolution par Critère d'Optimalité

Un schéma de résolution heuristique et simple de ce type du problème (3) est proposé par [BENDSOE-1995] et consite à mettre à jour, de manière itérative, les densités courantes \(\textbf{x}\) vers une nouvelle valeur \(\textbf{x}^{\textrm{new}}\) :

(4)\[\begin{split}x_e^{\textrm{new}} = \left\{ \begin{array}{lll} x_e^- & \textsf{si} & x_eB_e^{\eta} \le x_e^- \\ x_eB_e^{\eta} & \textsf{si} & x_e^- < x_eB_e^{\eta} < x_e^+ \\ x_e^+ & \textsf{si} & x_e^+ \le x_eB_e^{\eta} \\ \end{array} \right.\end{split}\]

où :

  • \(x_e^- = \max (0,x_e-m)\) est une borne inférieure pour la densité à l'itération courante

  • \(x_e^+ = \min (1,x_e+m)\) est une borne supérieure pour la densité à l'itération courante

  • \(m\) est la limite d'incrément de densité à chaque itération pour stabiliser la convergence

  • \(\eta\) est un coefficient d'amortissement (généralement \(\eta=0,5\))

Le terme \(B_e\) guidant la mise à jour de \(x_e\) est obtenu par la condition d'optimalité :

(5)\[B_e = \frac{-\dfrac{\partial f}{\partial x_e}}{\mathcal{L} \dfrac{\partial g}{\partial x_e}}\]

  • \(\dfrac{\partial f}{\partial x_e}\) est la sensibilité de la fonction objectif \(f\)

  • \(\dfrac{\partial g}{\partial x_e}\) est la sensibilité de la fonction contrainte \(g\)

  • \(\mathcal{L}\) est le multiplicateur de Lagrange à déterminer pour satisfaire la contrainte \(g\)

En dérivant les expressions des fonctions, la sensibilité de la fonction objectif (compliance), en supposant l'absence de forces dépendantes de la densité \((\frac{\partial \textbf{F}}{\partial x_e}=0)\), s'écrit :

(6)\[\frac{\partial f}{\partial x_e} = -p(x_e)^{p-1}(E_0-E_{\textrm{min}}) \textbf{u}_e^T.\mathbfcal{k}_0.\textbf{u}_e\]

La sensibilité de la fonction contrainte (fraction volumique) s'écrit :

(7)\[\frac{\partial g}{\partial x_e} = \frac{V_e}{V_0}\]

La difficulté étant alors de trouver la valeur de \(\mathcal{L}\) qui satisfait la contrainte \(g\). Celle-ci peut être calculée par un algorithme de dichotomie en initialisant des bornes inférieure \(\mathcal{L}^-\) et supérieure \(\mathcal{L}^+\) puis en choisissant la valeur milieu de l'intervalle. Une évaluation de la fonction contrainte \(g\) est alors faite et le processus est répété dans le demi intervalle ad hoc :

Initialisation des bornes

\(\mathcal{L}^- =0 \quad \mathcal{L}^+ =10^9\)

Tant que \(\frac{\mathcal{L}^+ - \mathcal{L}^-}{\mathcal{L}^+ + \mathcal{L}^-} > 10^{-3}\) :

\[\begin{split}\begin{array}{ll} \mathcal{L} & = (\mathcal{L}^- + \mathcal{L}^+)/2 \\ \textbf{B} & = \textsf{calculer selon (5)} \\ \textbf{x}^{\textrm{new}} & = \textsf{calculer selon (4)} \\ g(\textbf{x}^{\textrm{new}}) & = \textsf{calculer selon (3)} \\ \textsf{si } g(\textbf{x}^{\textrm{new}}) & > 0 & \\ \quad \mathcal{L}^- & = \mathcal{L} \\ \textsf{sinon} & \\ \quad \mathcal{L}^+ & = \mathcal{L} \\ \textsf{finsi} &\\ \end{array} \\\end{split}\]

Fin

À l'issue de la dichotomie on obtient la valeur de \(\mathcal{L}\) qui satisfait la contrainte sur la fraction volumique ainsi que la nouvelle topologie \(\textbf{x}^{\textrm{new}}\).

Filtrage de la sensibilité

Afin d'éviter l'effet de damier et diminuer la sensibilité des solutions au maillage, on applique une procédure de filtrage (ou lissage) sur le champ de desnité ou sur les champs de sensibilité.

Le filtrage de la sensibilité consiste à modifier la sensibilité \(\partial f / \partial x_e\) d'un élément \(e\) par une valeur moyenne pondérée des sensibilités des éléments voisins :

(8)\[\dfrac{\widehat{\partial f}}{\partial x_e} = \frac{1}{\textrm{max}(10^{-3},x_e)}\frac{\sum_{i \in N_e}\hat{H}_{ei}x_i\frac{\partial f}{\partial x_i}}{\sum_{i \in N_e}\hat{H}_{ei}}\]

Les coefficients \(\hat{H}_{ei}\) valent :

\[\hat{H}_{ei} = \left( 1 - \frac{\textrm{dist}(e,i)}{r_{\textrm{min}}} \right)^q V_i\]

et n'est définit que pour les \(N_e\) éléments \(i\) tels que \(\textrm{dist}(e,i) \le r_{\textrm{min}}\). Détaillons les nouvelles variables :

  • \(\textrm{dist}(e,i)\) la distance entre les centres des éléments \(e\) et \(i\)

  • \(V_i\) le volume de l'élément \(i\) (ou bien une autre quantitié pour pondérer)

  • \(r_{\textrm{min}}\) le rayon du filtre, au dela duquel l'opérateur \(\hat{H}_{ei}\) est nul

  • \(q\) un coefficient

Dans Cast3M, ce filtrage est réalisé grâce aux opérateurs MFIL (pour calculer le terme au numérateur) et NFIL (pour normaliser par le terme au dénominateur).

Notons que dans les articles des codes 99 lignes [SIGMUND-2001] et 88 lignes [ANDREASSEN-AL-2011], l'opérateur de filtrage utilisé correspond au cas où \(q=1\) et où tous les éléments ont un volume \(V_i=1\).

Filtrage des densités

A faire ...

Seuillage des densités

A faire ...