1. Introduction à l'optimisation topologique
1.1. 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é \(\mathbb{k}\), ce qui peut s'écrire :
Explicitons certains termes :
Par répartition optimale de la matière on entend la distribution spatiale de rigidité \(\mathbb{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. C'est l'ensemble des fonctions \(\chi_i(\mathbb{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 maximiser la raideur de la structure, 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 reformulé :
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
\(\mathbb{K}\) est la matrice de rigidité globale
\(V(\textbf{x})\) est le volume de la topologie \(\textbf{x}\) et \(V_e\) le volume de l'élément \(e\)
\(f\) la fraction volumique imposée
\(V_0\) le volume de domaine de conception \(\Omega\)
Remarquons que la compliance \(\psi(\textbf{x}) = \textbf{F}^T.\textbf{U}\) correspond aussi au travail des forces extérieures.
1.2. 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 inustriels et celle mise en oeuvre dans Cast3M via la procédure TOPOPTIM. Le lecteur intéressé pourra consulter de nombreux ouvrages sur le sujet dont :
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é \(\mathbb{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 :
Le problème d'optimisation de la compliance devient finalement :
avec :
\(\textbf{u}_e\) et \(\textbf{f}_e\) les vecteurs déplacements et forces de l'élément \(e\)
\(\mathbb{k}_e\) la matrice de rigidité de l'élément \(e\)
\(\mathbb{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\))
1.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}}\) :
avec :
\(m\) est une limite d'incrément de densité sur l'itération pour stabiliser la convergence
\(\eta\) est un coefficient d'amortissement (généralement \(\eta=0,5\))
\(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)\) la borne supérieure sur \(x_e\) pour respecter l'inégalité \(x_e^{\textrm{new}} \le 1\)
Le terme \(B_e\) guidant la mise à jour de \(x_e\) est obtenu par la condition d'optimalité :
\(\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\)
\(\lambda\) 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 :
La sensibilité de la fonction contrainte (volume) s'écrit :
La difficulté étant alors de trouver la valeur de \(\lambda\) qui satisfait la contrainte. Étant donné que la fonction contrainte \(\chi\) a une décroissance monotone avec \(\lambda\), on peut utiliser une dichotomie en initialisant des bornes inférieure \(\lambda^-\) et supérieure \(\lambda^+\) 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
\(\lambda^- =0 \quad \lambda^+ =100 000\)
Tant que \((\lambda^+ - \lambda^-) > 0,0001\) :
Fin
À l'issue de la dichotomie on obtient la valeur de \(\lambda\) qui satisfait la contrainte sur le volume ainsi que la nouvelle topologie \(\textbf{x}^{\textrm{new}}\).
1.4. 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}}\) :
L'opérateur de convolution \(\hat{H}_f\) vaut :
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
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\).
1.5. Illustration sur un cas mécanique
Une mise en donnée de l'algorithme d'optimisation précédent est fournie en annexe et disponible sur le site Cast3M.
Il s'agit d'optimiser la poutre en flexion présentée plus haut
Pour l'optimisation, on choisit :
une fonction objectif : la compliance \(\psi(\textbf{x}) = \textbf{U}^T.\mathbb{K}(\textbf{x}).\textbf{U}\)
une contrainte sur le volume : \(f=40\%\) du domaine de conception
les paramètres d'optimisation \(p=3\), \(\eta=0,5\), \(m=0,1\) et \(x_{\textrm{min}}=0,001\).
On initialise la topologie x avec des densités homogènes \(x_e=f\) afin de satisfaire la contrainte de volume.
Le volume cible est nommé vcib.
On calcule la matrice de filtrage kfil, intervenant dans l'équation (8) avec
l'opérateur MFIL.
Notons que pour cela, il est nécessaire de disposer du maillage mcg des centres de gravité du maillage ainsi que
du champ par points wg des volumes \(V_e\) de chaque éléments, exprimé sur ces centres de gravité.
Le champ des volumes élémentaires vole est obtenu grâce à l'opérateur INTG 'ELEM'
en intégrant un champ unitaire par élément.
Initialisation : topologie initiale et matrice de filtrage
44** Optimisation topologique
45* initialisation de la topologie
46x = MANU 'CHML' mod 'SCAL' fv 'GRAVITE' ;
47vini = MESU mail ;
48vcib = vini * fv ;
49* matrice de filtrage
50un = MANU 'CHML' mod 'SCAL' 1. 'GRAVITE' ;
51vole = INTG mod un 'ELEM' ;
52mcg = un POIN 'SUPERIEUR' 0. ;
53wg = PROI mcg vole ;
54kfil = MFIL wg rmin 1. 0. ;
On démarre ensuite une boucle d'optimisation limitée à 100 itérations.
On calcule alors rip, la matrice de rigidité pondérée \(\mathbb{K}(\textbf{x})\) de la
topologie courante selon la loi puissance de la méthode SIMP. Le comportement étant isotrope,
le module d'Young pénalisé yop de chaque élément vaut \(E_e=(x_e)^pE_0\) avec
\(E_0\) le module d'Young du matériau.
On résoud ensuite le problème mécanique \(\mathbb{K}(\textbf{x}).\textbf{U} =\textbf{F}\).
Pénalisation de la rigidité et résolution
55* boucle d'optimisation
56liso = PROG 0. 'PAS' 0.05 1. ;
57REPE b1 100 ;
58* pénalisation de la rigidité
59 yop = (x ** p) * yo ;
60 map = MATE mod 'YOUN' yop 'NU' nu ;
61* résolution du problème mécanique
62 rip = RIGI mod map ;
63 u = RESO (rip ET blo) f ;
On peut calculer la valeur psi de la fonction objectif \(\psi(\textbf{x}) = \textbf{F}^T.\textbf{U}\)
en remarquant que si celle-ci est égale au travail des forces extérieures, elle est donc aussi égale au travail des
efforts intérieurs et peut donc s'obtenir par :
où \(\sigma\) et \(\varepsilon\) désignent les champs de contraintes et déformations sig eps.
Le champ du double produit contracté \(\sigma:\varepsilon\) est obtenu grâce à l'opérateur
ENER et son intégrale par
INTG.
Le champ dpsi de sensibilité de la fonction objectif (6) s'exprime alors en fonction
de la matrice de Hooke \(\mathcal{C}_0\) du matériau plein :
Calcul de la fonction objectif et de sa sensibilités
67* fonction objectif : compliance = uT.K.u = Int(sig:eps)dV
68 eps = EPSI 'LINE' mod u ;
69 sig = ELAS mod map eps ;
70 psi = INTG mod (ENER mod sig eps) ;
71* sensibilité
72 sig0 = ELAS mod ma0 eps ;
73 ene0 = ENER mod eps sig0 ;
74 ene = CHAN ene0 mod 'GRAVITE' ;
75 dpsi = (-1. * p * (x ** (p - 1.)) * ene) ;
L'étape de filtrage de la sensibilité est réalisée en multipliant la matrice de filtrage kfil par le
champ par point xdpsi = x * dpsi représentant le produit \(x_f\dfrac{\partial \psi}{\partial x_f}\) dans
l'équation (8).
Filtrage de la sensibilité
76* filtrage de la sensibilité
77 xdpsi = x * dpsi ;
78 xdpsi0 = PROI mcg xdpsi ;
79 xdpsi1 = kfil * xdpsi0 ;
80 xdpsi = MANU 'CHML' mod 'REPA' 'SCAL' (EXTR xdpsi1 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
81 dpsi = xdpsi / x ;
La mise à jour de la topologie (passage du champ x à xnew) suivant le schéma (4)
est réalisée en suivant l'algorithme de dichotomie pour la recherche du
multiplicateur de Lagrange lmid qui nécessite une nouvelle boucle (limitée à 50 itérations).
La limitation d'incrément \(m\) et le recpect des bornes \(x_\textrm{min} \le x_e \le 1\) sont réaliséespace grâce aux opérateurs BORN et MASQ, ce dernier produisant un champ unitiare aux éléments respectant une inégalité et nul ailleur.
La vérification de la contrainte de volume est faite en calculant le volume vnew de chaque toplogie xnew
et en le comparant au volume cible vcib.
Optimisation (critère d'optimalité)
82* optimisation d'une nouvelle topologie (méthode du critère d'optimalité)
83 l1 = 0. ;
84 l2 = MAXI 'ABS' dpsi ;
85 REPE b2 50 ;
86 SI ((l2 - l1) < 1.E-4) ;
87 QUIT b2 ;
88 FINSI ;
89 lmid = 0.5 * (l1 + l2) ;
90 b = -1. * dpsi / lmid ;
91 xinf = BORN (x - m) 'MINIMUM' xmin ;
92 xsup = BORN (x + m) 'MAXIMUM' xmax ;
93 xnew = x * (b ** eta) ;
94 minf = (xnew - xinf) MASQ 'INFERIEUR' 0. ;
95 mmil = ((xnew - xinf) MASQ 'SUPERIEUR' 0.) * ((xnew - xsup) MASQ 'INFERIEUR' 0.) ;
96 msup = (xnew - xsup) MASQ 'SUPERIEUR' 0. ;
97 xnew = (xinf * minf) + (xnew * mmil) + (xsup * msup) ;
98 vnew = INTG mod xnew ;
99 SI (vnew > vcib) ;
100 l1 = lmid ;
101 SINON ;
102 l2 = lmid ;
103 FINSI ;
104 FIN b2 ;
Un affichage bilan de l'itération est fait, puis un cirtère d'arrêt de la boucle d'optimisation est proposé lorsque l'incrément maximal de densité est inférieur à 0,01
Fin de boucle et critère d'arrêt
105* bilan de l'itération
106 change = MAXI 'ABS' (x - xnew) ;
107 MESS 'It.' ' ' &b1 ' Obj.' ' ' c ' Change' ' ' change ' Fvol.' ' ' (vnew / vini) ' Dicho.' ' ' &b2 ;
108 SI (change < 0.01) ;
109 QUIT b1 ;
110 FINSI ;
111* préparation de la nouvelle itération
112 x = xnew ;
113FIN b1 ;
Les résultats de cette optimisation sont présentés dans l'animation ci-dessous qui montre les topologies (champs par éléments de densités) obtenues au cours des itérations. La topologie finale est atteinte après 42 itérations.
Animation des topologies au cours de l'optimisation (déformée x 1000)