# Une méthode de champ de phase pour l'optimisation topologique

Blaise Bourdin, Antonin Chambolle

<figure><img src="image_1.png" /></figure>

To cite this version:

Blaise Bourdin, Antonin Chambolle. Une méthode de champ de phase pour l'optimisation topologique. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812888

HAL Id: hal-01812888

https://hal.science/hal-01812888v1

Submitted on 11 Jun 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Une méthode de champ de phase pour l'optimisation topologique

Blaise Bourdin* — Antonin Chambolle**

* Department of Mathematics
Louisiana State University
Baton Rouge, LA 70803, USA
bourdin@math.lsu.edu

** CMAP, Ecole Polytechnique
CNRS (UMR 7641)
91128 Palaiseau Cedex, France
antonin.chambolle@polytechnique.fr

<footnote>RÉSUMÉ. Nous décrivons ici, en la simplifiant, une méthode d'optimisation topologique détaillée dans [BOU 03]. L'idée est de représenter la forme par un champ de phase v, et de régulariser le problème en ajoutant au critère a minimiser une pénalisation du périmètre total de la forme, ce qui revient à ajouter à l'énergie totale une fonctionnelle quadratique en v, relativement élémentaire (bien que non convexe). Nous expliquons ensuite sommairement comment on peut résoudre le problème, et montrons quelques résultats numériques en dimension 2.</footnote>

<footnote>ABSTRACT. We describe in this note a simplified version of our method of topology optimization proposed in [BOU 03]. The idea is that the shape is represented by a phase-field v, and the optimization problem is regularized by a penalization of the total perimeter of the shape. This amounts to adding to the energy a quadratic functional of v, which is quite simple (although nonconvex). We explain briefly our approach to solve the problem, and show some numerical results in 2D.</footnote>

<footnote>MOTS-CLÉS : optimisation topologique, pénalisation du périmètre, interface diffuse, champ de phase, $\Gamma$-convergence.</footnote>

<footnote>KEYWORDS: topology optimization, perimeter penalization, diffuse interface, phase-field, $\Gamma$-convergence.</footnote>
# 1. Introduction

Le calcul de formes optimales est un problème difficile, à la fois parce qu'il n'est pas évident de trouver une représentation aisément manipulable d'une forme (d'autant moins si sa topologie est a priori inconnue), mais aussi parce qu'il s'agit presque toujours d'un problème mal posé : les phénomènes d'homogénéisation font que l'optimum est rarement une “forme” au sens classique. Le but des méthodes de “champ de phase” est de remédier simultanément à ces deux problèmes : d'un côté, en ajoutant au critère à minimiser une pénalisation du périmètre global de la forme (qui interdit l'homogénéisation), et de l'autre, en approchant le critère global ainsi obtenu par une énergie dépendant d'une “phase” $v \in [0, 1]$, qui va tendre vers la fonction caractéristique de la forme inconnue lorsque qu'un paramètre $\varepsilon$ tend vers zéro. L'idée de pénaliser le périmètre a été introduite, en optimisation de forme, dans [AMB 93] et [HAB 96], et notamment étudiée par J. Petersson [PET 99]. Pour un exposé général sur l'optimisation topologique nous renvoyons à [BEN 02].

Nous présentons ici de façon simplifiée cette méthode, appliquée dans [BOU 03] au problème d'optimisation de la compliance d'un matériau soumis à des forces volumiques et éventuellement de pression, sur une partie de la surface inconnue. Dans ce cas, le champ de phase doit se concentrer sur trois phases principales, le vide, l'objet, et un liquide fictif qui exerce une pression au contact avec l'objet. Pour simplifier, nous ne considérons ici qu'un système à deux phases (solide/vide). Récemment, la même méthode a été appliquée par Wang et Zhou à des problèmes du même type [WAN 03], tandis que Burger et Stainko viennent de montrer comment l'adapter au cas de la minimisation du volume lorsque l'on impose à la contrainte mécanique de rester locallement bornée, avec des résultats remarquables (le cas d'un critère de type Von Mises est notamment considéré) [BUR 04].

# 2. L'approximation du périmètre

La méthode est basée sur la remarque suivante, qui découle d'un théorème de Γ-convergence d'abord démontré par Modica et Mortola [MOD 77a, MOD 77b] (voir aussi [BAL 90, ALB 00]). Dans la suite, $\Omega \subset \mathbb{R}^d$ ($d = 2$ ou 3) désignera la “boîte” dans laquelle on cherche à calculer un domaine optimal. Soit $C(v)$ un critère continu dépendant de la variable $v(x) \in L^1(\Omega)$, et considérons les fonctionnelles

$$ \mathcal{P}_{\varepsilon}(v) = \frac{\varepsilon}{2} \int_{\Omega} |\nabla v(x)|^2 dx + \frac{1}{\varepsilon} \int_{\Omega} W(v(x)) dx, \text{ et} $$

$$ \mathcal{E}_{\varepsilon}(v) = C(v) + \lambda \mathcal{P}_{\varepsilon}(v), $$

où $\lambda > 0$ est un paramètre réel (petit), et $W(t)$ un potentiel s.c.i. qui vérifie $W(0) = W(1) = 0$, $W(t) > 0$ si $t \notin \{0, 1\}$, et $W(t) \ge c|t|^2 - c'$.

Alors, si pour tout $\varepsilon$ petit, $v_{\varepsilon}$ minimise $\mathcal{E}_{\varepsilon}$, il existe une sous-suite $(v_{\varepsilon_j})_{j \ge 1}$ et un ensemble $E \subset \Omega$ tel que $v_{\varepsilon_j}(x) \to \chi_E(x)$ (la caractéristique de $E$) p.p. dans $\Omega$, et
$E$ minimise $C(\chi_F) + \lambda c_W$ Périmètre($F$) parmi tous les ensembles $F$ de périmètre fini dans $\Omega$. De plus, si pour $\varepsilon > 0$, $v_\varepsilon$ est un minimiseur de $\mathcal{E}_\varepsilon$ sous contrainte $\int_\Omega v \, dx = V$, alors l'ensemble limite $E$ minimise $C(\chi_F) + \lambda c_W$ Périmètre($F$) sous la contrainte de volume $|F| = V$. La constante $c_W$ est donnée par $c_W = \int_0^1 \sqrt{2W(t)} \, dt$.

On voit qu'on peut donc trouver de bonnes approximation d'un minimiseur de $L(\chi_F)$, plus une pénalisation du périmètre, en minimisant $\mathcal{E}_\varepsilon$ pour $\varepsilon$ petit. Cette fonctionnelle a l'avantage de ne dépendre que d'une fonction à valeurs réelles. Un "bon" choix de $W$ est par exemple

$$W(t) = \begin{cases} t(t-1)/2 & \text{si } 0 \le t \le 1, \\ +\infty & \text{sinon,} \end{cases} \quad [1]$$

pour lequel on a $c_W = \pi/8$.

Ce résultat présente un inconvénient : il n'est pas garanti du tout que l'ensemble limite $E$ soit régulier, ni même ouvert : en pratique, cela va dépendre fortement du critère $C$. En fait, $E$ vit dans la classe des "ensembles de périmètre fini" (voir par exemple [GIU 84]) qui contient quelques ensembles assez pathologiques tels que $\bigcup_{i=1}^\infty B(x_i, 2^{-i})$ où $(x_i)_{i=1}^\infty$ est la suite des points à coordonnées rationnelles de $\Omega...$ On sait seulement dans certain cas (voir [AMB 93, CHA 03]) que l'ensemble minimisant est en réalité plus régulier. Dans la suite, nous noterons $\mathcal{P}(E)$ le périmètre généralisé de l'ensemble $E$ (que l'on peut définir comme la plus petite limite des périmètres des $E_n$ pour toutes les suites d'ensembles réguliers $E_n$ tendant vers $E$).

Le résultat ci-dessus doit être un peu précisé : si aucune condition au bord n'est imposée à $v_\varepsilon$, à la limite, $\mathcal{P}(E)$ ne pénalise que la partie du périmètre qui est dans $\Omega$ (grosso modo, la surface totale de $\partial E \cap \Omega$). Si par contre on impose sur une partie de la frontière $\Gamma \subset \partial \Omega$ une condition de Dirichlet $v_\varepsilon = 0$, alors, à la limite, $\mathcal{P}(E)$ pénalise aussi la surface de $\partial E \cap \Gamma$. Enfin, si une condition $v_\varepsilon = 1$ est imposée sur $\Gamma$, $\mathcal{P}(E)$ pénalise $\Gamma \setminus \partial E$.

### 3. Le problème faible

Nous considérons le problème de minimisation de la compliance d'une structure linéaire-élastique, sous une contrainte de volume. On se donne $\Gamma^0, \Gamma^f \in \partial \Omega$ deux parties de la frontière à distance strictement positive l'une de l'autre, et, sur $\Gamma^f$, une force $f(x)$. Etant donné $E \subset \Omega$, on définit la compliance de $E$ par

$$C(E) = - \inf_{\mathbf{u}=0 \text{ sur } \Gamma^0} \int_E \mathbf{A} e(\mathbf{u}) : e(\mathbf{u}) \, dx - 2 \int_{\Gamma^f} \mathbf{f} \cdot \mathbf{u} \, d\sigma.$$

où $\mathbf{A}$ est la loi de Hooke d'un matériau linéaire-élastique, $\mathbf{u}$ est un déplacement (de $\Omega$ à valeurs dans $\mathbb{R}^d$), et $e(\mathbf{u})$ le gradient symétrisé de $\mathbf{u}$. Si $E$ est un ensemble régulier, $C$ représente la compliance (le travail de la force $\mathbf{f}$ à l'équilibre) de la structure de configuration de référence $E$. En particulier, $C(E)$ est infinie lorsque $\Gamma^f \not\subset \partial E$. L'avantage de cette formulation est qu'elle est correctement définie, y compris
lorsque E est un simple ensemble de périmètre fini. On peut donc chercher un minimum de $C(E) + \lambda c_W P(E)$ dans cette classe, sous contrainte de volume $|E| \le V$ ($V \in ]0, |\Omega|]$), et on montre aisément que ce minimum existe. Dans [BOU 03], nous traitons aussi le cas de forces internes, et de forces de pression s'exerçant sur une partie a priori inconnue de la frontière de E.

Cependant, la fonctionnelle $\mathcal{C}$ que nous venons de définir, vue comme une fonction de l'ensemble caractéristique $\chi_E$, n'est que s.c.i. dans $L^1(\Omega)$ (et non continue). On ne peut donc approcher directement le problème de minimisation de $\mathcal{C}(E) + \lambda c_W P(E)$ par des problèmes de la forme $\mathcal{C}(v) + \lambda P_\varepsilon(v)$, où $\mathcal{C}$ est étendue (de façon évidente) à tout $L^1(\Omega)$. On est obligé d'introduire dans $\Omega \setminus E$ un matériau fictif, de loi de Hooke (par exemple) $\delta \mathbf{A}$, où $\delta > 0$ est un petit paramètre. On définit alors

$$ \mathcal{C}_{\delta}(E) = \inf_{\mathbf{u}=0 \text{ sur } \Gamma^0} \int_{\Omega} ((1-\delta)\chi_E + \delta)\mathbf{A}e(\mathbf{u}) : e(\mathbf{u}) \, dx - 2 \int_{\Gamma^f} \mathbf{f} \cdot \mathbf{u} \, d\sigma. $$

La nouvelle fonctionnelle $\mathcal{C}_\delta$ est toujours finie et il existe toujours un déplacement minimisant $\mathbf{u} \in H^1(\Omega; \mathbb{R}^d)$. Toujours grâce à un résultat de $\Gamma$-convergence [BOU 03], on montre que si $E_\delta$ minimise $\mathcal{C}_\delta(E) + \lambda c_W P(E)$ pour $\delta$ petit, alors on peut extraire une sous-suite $(E_{\delta_j})_{j \ge 1}$ et trouver un ensemble de périmètre fini $E$ tel que $\chi_{E_{\delta_j}} \to \chi_E$ dans $L^1(\Omega)$ quand $j \to \infty$: dans ce cas, $E$ est alors forcément un minimiseur de $\mathcal{C}(E) + \lambda c_W P(E)$. On voit qu'on a introduit une "bonne" approximation du problème faible.

L'avantage de cette nouvelle formulation est qu'il est évident de prolonger $\mathcal{C}_\delta$ à tout $L^1(\Omega)$, en une fonctionnelle continue. On pose, pour $v \in L^1(\Omega)$,

$$ \mathcal{C}_{\delta}(v) = \inf_{\mathbf{u}=0 \text{ sur } \Gamma^0} \int_{\Omega} ((1-\delta)S(v) + \delta)\mathbf{A}e(\mathbf{u}) : e(\mathbf{u}) \, dx - 2 \int_{\Gamma^f} \mathbf{f} \cdot \mathbf{u} \, d\sigma. $$

où $S(t)$ est une fonction régulière, croissante, telle que $S(t) = 0$ si $t \le 0$ et $S(t) = 1$ si $t \ge 1$. La continuité de $\mathcal{C}_\delta$ n'est pas difficile à démontrer. On peut alors approcher l'énergie $\mathcal{C}_\delta(E) + \lambda c_W P(E)$ par $\mathcal{C}_\delta(v) + \lambda P_\varepsilon(v)$, $\varepsilon > 0$. On a alors, si $v_{\varepsilon,\delta}$ minimise cette fonctionnelle, un double processus de convergence : si $\varepsilon \to 0$, $v_{\varepsilon,\delta} \to \chi_{E_\delta}$ (à une sous-suite près), et $E_\delta$ minimise $\mathcal{C}_\delta(E) + \lambda c_W P(E)$ ; donc si $\delta \to 0$, $E_\delta \to E$ (à une sous-suite près), où $E$ minimise le problème faible $\mathcal{C}(E) + \lambda c_W P(E)$.

## 4. L'implémentation numérique

Pour minimiser $\mathcal{C}_\delta(v) + \lambda P_\varepsilon(v)$ par rapport à $v$, sous la contrainte de volume $\int_\Omega v \, dx = V$, on discrétise le problème sur une base d'éléments finis, puis on implémente un schéma de descente de gradient, semi-implicite en temps. La dérivée de $\mathcal{C}_\delta$ est calculée explicitement : en effectuant des petites variations, on trouve

$$ \langle D\mathcal{C}_\delta(v), w \rangle = \lim_{t \to 0} \frac{\mathcal{C}_\delta(v + tw) - \mathcal{C}_\delta(v)}{t} = - \int_\Omega (1 - \delta)S'(v)\mathbf{A}e(\mathbf{u}) : e(\mathbf{u})w \, dx \quad [2] $$
où **u** est la position d'équilibre pour **v** (ce calcul peut être un peu plus compliqué pour d'autres critères et nécessiter le calcul d'un état adjoint).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Le piston</figcaption></figure>

On résout alors (dans le cas où $W$ est l'exemple donné en [1]), d'où $W'(t) = 1/2 - t$ si $|t| < 1$,

$$ \frac{v^{n+1} - v^n}{\delta t} = \lambda \epsilon \Delta v^{n+1} - \frac{\lambda}{\epsilon} \left( \frac{1}{2} - v^{n+1} \right) - D C_{\delta}(v^n) + \theta \quad [3] $$

qui est bien posé lorsque $\delta t \le \epsilon/\lambda$. Le paramètre $\theta$ sert à contrôler la contrainte $\int_{\Omega} v^{n+1} dx = V$ et doit être ajusté à chaque itération. Les valeurs de $v^{n+1}$ sont ensuite ramenées entre 0 et 1 par troncature (cette étape de reprojection est dictée par le choix particulier d'un $W$ infini hors de $[0, 1]$, un autre choix, tel que $W(t) = t^2(1 - t^2)/4$, conduirait à une évolution légèrement différente de [3], mais sans reprojection).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Le pont</figcaption></figure>

On itère, en partant initialement d'une valeur de $\epsilon$ relativement grande, et en la réduisant peu à peu. On peut optimiser le maillage au fur et à mesure, en le raffinant là où le gradient de $v$ est le plus élevé. La valeur de $\epsilon > 0$ est liée à l'épaisseur de la zone de transition (si $W$ est la fonction [1], on montre que l'interface diffuse est d'épaisseur $\simeq \epsilon\pi$), et il est important que la taille du maillage dans cette zone soit inférieure.

Nous montrons deux exemples de résultats calculés avec cette méthode : la figure 1, extraite de [BOU 03], représente un piston, fixé au petit rectangle noir, et
subissant une pression sur sa partie inférieure. Cet exemple a été calculé avec un
code d'élément finis développé par Blaise Bourdin. La figure 2 représente un pont,
en configurations de référence et déformée. L'objet, posé aux deux extrémité, est sou-
mis à une force au milieu du bord inférieur (la moitié du pont seulement a été cal-
culée, puis symétrisée). Cet exemple a été calculé sous FreeFem++, un logiciel gra-
tuit d'éléments finis 2D disponible sur www.freefem.org. Le programme, court et
relativement simple, a été écrit avec la collaboration de deux étudiants de l'Ecole Po-
lytechnique, Arnaud Anantharaman et Alain Griveau. Il implémenta exactement les
équations [2] et [3] (W est donnée par [1], et on a choisi S(t) = t² si t ∈ [0, 1]).

**5. Bibliographie**

[ALB 00] ALBERTI G., « Variational models for phase transitions, an approach via Γ-convergence », *Calculus of variations and partial differential equations* (Pisa, 1996), p. 95–114, Springer, Berlin, 2000.

[AMB 93] AMBROSIO L., BUTTAZZO G., « An optimal design problem with perimeter penalization », *Calc. Var. Partial Differential Equations*, vol. 1, n° 1, 1993, p. 55–69.

[BAL 90] BALDO S., « Minimal interface criterion for phase transitions in mixtures of Cahn-Hilliard fluids », *Ann. Inst. H. Poincaré Anal. Non Linéaire*, vol. 7, n° 2, 1990, p. 67–90.

[BEN 02] BENDSØE M. P., SIGMUND O., *Topology Optimization*, Springer, Berlin, 2002.

[BOU 03] BOURDIN B., CHAMBOLLE A., « Design-dependent loads in topology optimization », *ESAIM Control Optim. Calc. Var.*, vol. 9, 2003, p. 19–48 (electronic).

[BUR 04] BURGER M., STAINKO R., « Phase-field relaxation of topology optimization with local stress constraints », rapport n° 2004-35, 2004, J. Kepler University Linz, SFB F013.

[CHA 03] CHAMBOLLE A., LARSEN C. J., « $C^\infty$ regularity of the free boundary for a two-dimensional optimal compliance problem », *Calc. Var. Partial Differential Equations*, vol. 18, n° 1, 2003, p. 77–94.

[GIU 84] GIUSTI E., *Minimal surfaces and functions of bounded variation*, Birkhäuser, Boston, 1984.

[HAB 96] HABER R., JOG C., BENDSØE M., « A new approach to variable-topology shape design using a constraint on the perimeter », *Struct. Multidisc. Optim.*, vol. 11, 1996, p. 1–12.

[MOD 77a] MODICA L., MORTOLA S., « Il limite nella Γ-convergenza di una famiglia di funzionali ellittici », *Boll. Un. Mat. Ital. A (5)*, vol. 14, n° 3, 1977, p. 526–529.

[MOD 77b] MODICA L., MORTOLA S., « Un esempio di Γ⁻-convergenza », *Boll. Un. Mat. Ital. B (5)*, vol. 14, n° 1, 1977, p. 285–299.

[PET 99] PETERSSON J., « Some convergence results in perimeter-controlled topology optimization », *Comput. Methods Appl. Mech. Engrg.*, vol. 171, n° 1-2, 1999, p. 123–140.

[WAN 03] WANG M. Y., ZHOU S. W., « Phase transition : a variational method for structural topology optimization », Preprint Chinese University of HongKong, 2003.