Algorithme numérique d'optimisation de formes

Plus encore que le théorème d'existence de solutions relaxées, l'intérêt de cette méthode d'homogénéisation provient de l'algorithme numérique d'optimisation qu'on en déduit. Le principe de la méthode proposée est de calculer une forme composite optimale pour la formulation homogénéisé (6), plutôt que d'essayer de trouver une forme classique quasi-optimale pour la formulation originale (4). En effet, rappelons que ce problème (4) est mal posé et n'admet pas en général de solution classique. Numériquement, cela veut dire que la minimisation de (4) est très difficile à cause de la présence de nombreux minima locaux : les résultats numériques dépendent alors fortement des conditions initiales et du maillage choisis. Au contraire le problème homogénéisé (6) admet un minimum global. De plus, l'homogénéisation transforme un difficile problème d'optimisation discrète (du type 0/1 selon qu'au point $x\in\Omega$ il y a, ou non, du matériau) en un problème beaucoup plus simple d'optimisation continue (la densité de matériau varie continuement entre 0 et 1). Les calculs sont donc effectués sur un maillage fixe du domaine de travail $\Omega$ sur lequel la forme optimale sera capturée par notre algorithme.

Le problème homogénéisé (6) d'optimisation de formes est, comme on l'a vu, une double minimisation par rapport aux paramètres de formes $(\theta,A^*)$ et aux tenseurs des contraintes statiquement admissibles $\sigma$. Notre algorithme, dit des ``directions alternées" (cf. [3]), se propose de minimiser alternativement et itérativement dans chacune de ces variables :

• intialisation de la forme $(\theta_0, A_0^*)$
• itérations $n\geq1$ jusqu'à convergence
  • étant donnée une forme $(\theta_{n-1}, A_{n-1}^*)$, on calcule les contraintes $\sigma_n$ par résolution d'un problème d'élasticité linéaire (par une méthode d'éléments finis)
  • étant donné ce tenseur des contraintes $\sigma_n$, on calcule les nouveaux paramètres de forme $(\theta_{n}, A_{n}^*)$ avec les formules explicites d'optimalité (8)-(11) faisant intervenir $\sigma_n$.

Figure 2: Chargement d'une poutre 2-D.

Comme chaque étape de ce procédé itératif est une minimisation partielle, la valeur de la fonction objectif décroit toujours, ce qui assure la convergence de la méthode. La partie d'optimisation des paramètres de formes est locale et facile grâce aux formules explicites. Tout l'effort de calcul porte donc sur la résolution d'une succession de problèmes d'élasticité linéaire. De nombreux calculs ont été effectués avec cet algorithme qui converge très vite indépendemment du choix initial $(\theta_0, A_0^*)$ et qui est stable par raffinement du maillage. Nous présentons ici un exemple du type "poutre en flexion" (cf. figure 2).

Figure 3: Poutre composite optimale.

Evidemment, notre algorithme calcule des formes composites (dont on trace la densité sur la figure 3) alors que dans la pratique on préfère obtenir de vraies formes classiques. Le problème est donc de faire disparaitre les zones composites (en gris sur les figures) pour retrouver une forme nette (image blanc et noir). Pour cela on introduit une technique de pénalisation des densités intermédiaires qui force la densité à ne prendre que les valeurs 0 ou 1. La stratégie est la suivante : après convergence de l'algorithme vers une forme composite optimale, on effectue encore quelques itérations de l'algorithme en pénalisant les densités intermédiaires. On obtient alors une forme classique capturée sur le maillage fixe. La raison du succès (un peu surprenant au premier abord, cf. figure 4) de cette étape de pénalisation est que la forme composite optimale est aussi caractérisée par une microstructure de perforations à l'échelle sous-maille (c'est-à-dire homogénéisée). Le fait de pénaliser les zones composites conduit l'algorithme à reproduire, au moins partiellement, cette microstructure à l'échelle du maillage. La pénalisation apparait donc comme une projection des formes composites sur des formes classiques. Bien sûr, cette forme classique quasi-optimale ainsi obtenue par projection est très dépendante du maillage. Plus ce dernier sera fin et plus elle incluera de détails liés à l'échelle du maillage. Mais du moins, cette forme classique reste proche du minimum global atteint par la forme composite optimale.

Une des propriétés les plus intéressantes de cette méthode est qu'elle permet d'obtenir des topologies très complexes de formes optimales sans aucune connaissance a priori de la solution. Notre algorithme permet donc de capturer sur un maillage fixe une forme optimale dont la topologie n'est pas fixée ni dictée par le choix des conditions initiales. Bien sûr, les contours de la forme ne sont pas très précis car ils ne coincident pas avec le maillage. Dans ces conditions cette méthode d'optimisation topologique apparait comme un préprocesseur qui peut s'utiliser avant un algorithme classique d'optimisation par variations de frontières. Il n'y a donc pas concurrence, mais plutôt complémentarité, entre ces deux approches.

Figure 4: Poutre optimale après pénalisation.

Nous renvoyons à [3], [4] pour tous les (nombreux) détails techniques numériques (par exemple, il faut tronquer la densité près de 0 pour éviter une singularité de la loi de Hooke homogénéisée). On peut aussi réintroduire les contraintes technologiques de faisabilité dans l'étape de pénalisation (voir par exemple [7], [29]). Par rapport au travail initial de Bendsoe et Kikuchi [8] qui utilisait comme microstructure une cellule périodique carrée percé d'un trou rectangulaire, notre approche se distingue (entre autres) par l'utilisation de microstructures optimales, à savoir les laminés séquentiels. Il existe de nombreuses généralisations du problème modèle présenté ici. Citons entre autres l'optimisation multi-chargements (avec plusieurs équations d'état), l'optimisation de fréquences propres de vibration, l'optimisation de modèles de plaques (voir par exemple les références dans [7]).

Figure 5: Pylône optimal en 3-D (déplacements nuls aux quatre coins en bas, forces ponctuelles dirigées vers le bas aux points d'ancrage des cables).

Remerciements Le travail décrit ici a été effectué en grandes parties avec la collaboration de E. Bonnetier, G. Francfort, F. Jouve, et R. Kohn. Qu'ils soient donc remerciés pour leur aide et leur amitié. Je remercie particulièrement François Jouve qui a développé le code avec lequel les résultats numériques présentés ici ont été obtenus.