Éléments théoriques de base

Généralités

L'optimisation topologique consiste à chercher 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 & \psi(\mathbfcal{k}) & \\ \textsf{tel que} \quad & \chi_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 une fonction \(\psi\) 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 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 \(\chi_i(\mathbfcal{k}) = 0\). En effet, afin d'éviter les solutions triviales, il est nécessaire d'imposer des contraintes. Par exemple, on peut chercher à minimiser le volume d'une pièce sous contrainte que celle-ci ne se déplace par trop (sans quoi la solution à volume nul sera optimale).

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

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

Plusieurs choix de fonctions \(\psi\) et \(\chi_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 de problème d'optimisation en mécanique des structures 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. En introduisant une discrétisation spatiale (maillage) et une variable de conception discrète \(\textbf{x}\), le problème d'optimisation topologique peut être plus clairement formulé :

(2)\[\begin{split}\min_{\textbf{x}} \quad & \psi(\textbf{x}) = \textbf{F}^T.\textbf{U}(\textbf{x}) = \textbf{U}^T(\textbf{x}).\mathbfcal{K}(\textbf{x}).\textbf{U}(\textbf{x}) \\ \textsf{tel que} \quad & \chi(\textbf{x}) = V(\textbf{x}) - fV_0 = \sum_{e=1}^N x_eV_e - fV_0 = 0 \\ & \mathbfcal{K}(\textbf{x}).\textbf{U}(\textbf{x}) =\textbf{F} \\ & x_e \in \{0;1\} \\\end{split}\]

Ce problème de minimiser de la compliance est équivalent à maximiser la raideur de la structure. Détaillons les variables :

• \(\textbf{x}\) est le vecteur des variables de conception \(x_e\), traduisant la présence de matière \((x_e=1)\) ou bien l'absence de matière \((x_e=0)\) dans l'élément \(e\)

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

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

• \(V(\textbf{x})\) est le volume de la topologie \(\textbf{x}\) et \(V_e\) le volume de l'élément \(e\)

• \(f\) est la fraction volumique imposée

• \(V_0\) est le volume du domaine de conception \(\Omega\)

La méthode SIMP

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-SIGMUND-2004] qui détaillent rigouresement la théorie derrière l'optimisation topologique.

• L'article pédagogique [SIGMUND-2001] qui présente une implémentation sur Matlab en 99 lignes d'un algorithme d'optimisation topologique. La procédure TOPOPTIM de Cast3M, ainsi que l'exemple utilisé dans ce document en sont grandement inspirés.

Les principales idées de la méthode sont les suivantes :

• Introduire des variables de conception continues \(x_e \in [0;1]\), appelées aussi densités.

• Pénaliser la rigidité \(\mathbfcal{K}\) en fonction de \(\textbf{x}\) par une loi puissance afin d'éviter la présence de densités intermédiaires. La matrice de rigidité de l'élément \(e\) vaut ainsi :

\[\mathbfcal{k}_e=(x_e)^p\mathbfcal{k}_0\]

Le problème d'optimisation de la compliance :eq:`opti_topo_2` devient finalement :

(3)\[\begin{split}\min_{\textbf{x}} \quad & \psi(\textbf{x}) = \textbf{U}^T(\textbf{x}).\mathbfcal{K}(\textbf{x}).\textbf{U}(\textbf{x}) = \sum_{e=1}^N (x_e)^p \quad \textbf{u}_e^T.\mathbfcal{k}_0.\textbf{u}_e \\ \textsf{tel que} \quad & \chi(\textbf{x}) = \sum_{e=1}^{N}x_eV_e - fV_0 = 0 \\ & \mathbfcal{K}(\textbf{x}).\textbf{U}(\textbf{x}) =\textbf{F} \\ & 0 < x_{\textrm{min}} \le x_e \le 1 \\\end{split}\]

avec :

• \(\textbf{u}_e\) et \(\textbf{f}_e\) les vecteurs déplacements et forces de l'élément \(e\)

• \(\mathbfcal{k}_e\) la matrice de rigidité de l'élément \(e\)

• \(\mathbfcal{k}_0\) la matrice de rigidité du matériau plein

• \(x_{\textrm{min}}\) une densité minimale non nulle (pour éviter les singularités)

• \(p\) le paramètre de pénalisation (en général \(p=3\))

Résolution du problème 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}\]

avec :

• \(x_e^- = \max (x_{\textrm{min}},x_e-m)\) une borne inférieure pour respecter l'inégalité \(x_{\textrm{min}} \le x_e^{\textrm{new}}\)

• \(x_e^+ = \min (1,x_e+m)\) une borne supérieure pour respecter l'inégalité \(x_e^{\textrm{new}} \le 1\)

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

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

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

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

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

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

• \(\mathcal{L}\) est un multiplicateur de Lagrange pour satisfaire la contrainte de volume \(\chi\)

En dérivant les expressions des fonctions, la sensibilité de la fonction objectif (compliance), en l'absence de forces dépendantes de la densité, s'écrit :

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

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

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

La difficulté étant alors de trouver la valeur de \(\mathcal{L}\) qui satisfait la contrainte. Étant donné que la fonction contrainte \(\chi\) a une décroissance monotone avec \(\mathcal{L}\), on peut utiliser une 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 \(\chi\) 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}^+ =100000000\)

Tant que \((\mathcal{L}^+ - \mathcal{L}^-) > 0,0001\) :

\[\begin{split}\begin{array}{ll} \mathcal{L} & = (\mathcal{L}^- + \mathcal{L}^+)/2 \\ \textbf{x}^{\textrm{new}} & = \textsf{actualiser } \textbf{x} \textsf{ selon (4)} \\ \textsf{si } \chi(\textbf{x}) & > 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 le volume 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) du champ de sensibilité. Sur chaque élément \(e\) la sensibilité de la compliance est remplacée par une valeur moyenne pondérée des sensibilités calculées sur les éléments voisins \(f\) dans un rayon \(r_{\textrm{min}}\) :

(8)\[\dfrac{\widehat{\partial \psi}}{\partial x_e} = \frac{1}{x_e}\dfrac{1}{\sum_{f=1}^{N_e}\hat{H}_f}\sum_{f=1}^{N_e}\hat{H}_fx_f\frac{\partial \psi}{\partial x_f}\]

L'opérateur de convolution \(\hat{H}_f\) vaut :

\[\hat{H}_f = \left( 1 - \frac{\textrm{dist}(e,f)}{r_{\textrm{min}}} \right)^q V_f\]

et n'est définit que pour les \(N_e\) éléments \(f\) tels que \(\textrm{dist}(e,f) \le r_{\textrm{min}}\), avec :

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

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

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

• \(q\) un coefficient

Dans Cast3M, ce filtrage est réalisé grâce à l'opérateur MFIL.

Notons que dans l'article 99 lignes de [SIGMUND-2001] l'opérateur de filtrage utilisé correspond au cas où \(q=1\) et où tous les éléments ont un volume \(V_f=1\).