# Optimisation topologique de problèmes fluidiques réactifs

Yann Favennec, Florian Dugast, Christophe Josset

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

## To cite this version:

Yann Favennec, Florian Dugast, Christophe Josset. Optimisation topologique de problèmes fluidiques réactifs. 24ème Congrès Français de Mécanique CFM 2019, Aug 2019, Brest, France. hal-02870118

HAL Id: hal-02870118

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

Submitted on 26 Jul 2021

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 CC BY 4.0 - Attribution - International License
# Optimisation topologique de problèmes fluidiques réactifs

Y. FAVENNEC<sup>a,*</sup>, F. DUGAST<sup>a</sup>, C. JOSSET<sup>a</sup>

a. Laboratoire de Thermique et Energie de Nantes, Polytech Nantes, Université de Nantes
yann.favennec@univ-nantes.fr

...

## Résumé :

L'optimisation topologique des transferts de chaleur et de masse et des échangeurs réactifs consiste à rechercher une géométrie optimale, au sens d'une distribution discrète solide/fluide, afin de satisfaire, au mieux, un objectif physique. Mathématiquement, il s'agit d'un problème de minimisation sous contraintes. Nous proposons un algorithme de descente basé sur une courbe de niveaux (level-set). Le gradient de la fonction de coût (qui représente l'objectif) est dérivé par états adjoints. Les modèles fluïdique et réactif sont résolus par des schémas de lattice Boltzmann, en version MRT (multi-relaxation time model) pour le problème convectif, et en version SRT (single-relaxation time model) pour la partie réactive. Quelques cas bi-dimensionnels d'application sont présentés dans ce congrès, autant pour la partie transfert de chaleur (par exemple maximisation de chaleur évacuée, ou minimisation d'une température moyenne) que pour la partie échangeur réactif (maximisation d'une réaction chimique dans un réacteur). Des contraintes de porosité ou de perte de charge sont également adjointes. On verra notamment que les géométries optimales obtenues n'auraient pas pu être pressenties. On s'appuiera sur une extension des développements présentés dans [4].

## Abstract :

Topology optimization of heat and mass transfer and of reactive devices consists in searching for an optimal geometry, in the sense of a discrete solid/fluid distribution, in order to satisfy, as most as possible, a given physical objective. Mathematically, such an objective is formulated as a (constrained) minimization problem. We here propose a gradient-based algorithm based on a level-set function. The cost function gradient is derived thanks to the adjoint-state approach. Both the fluid-flow model and the reaction model are based on the lattice Boltzmann method, with its multi-relaxation time model (MRT) version for the fluid flow, and the single-relaxation time model for the concentration model. A couple of two-dimensional applications are dealt with, on both heat transfer applications (for example maximization of some evacuated heat, or minimization of a mean temperature) and on reaction ones. Constraints on porosity and pressure drop are also integrated. The found optimal geometries could not have been predicted. These results constitute an extension of [4].
# 1 Introduction et résolution du problème d'optimisation

En optimisation topologique de problèmes fluidiques, il s'agit de trouver la meilleure distribution – ou du moins une distribution suffisamment bonne – pour qu'un certain objectif (un critère) soit satisfait. La distribution est représentée par une variable de design, et l'objectif est défini en terme de fonction de coût $\mathcal{J}$. La topologie $\mathcal{T}$, dont la construction est explicitée plus bas, est mise à jour à partir de l'évolution d'une courbe de niveau $\Psi$, elle-même basée sur le gradient spatial de la fonction de coût. Les contraintes d'inégalité (portant par exemple sur les pertes de charges) sont intégrées à la fonction de coût, par le biais de fonctions continues et dérivables, à valeurs élevées à l'activation des contraintes, et faibles sinon ; des fonctions exponentielles sont typiquement utilisées. On cherche :

$$ \mathcal{T}^* = \arg \min_{\mathcal{T}} \mathcal{J}(\cdot) \qquad (1) $$

Le domaine d'étude est divisé en $N$ éléments. Une variable de design signée est assignée à chaque élément, $\alpha(\boldsymbol{x}) \in \mathbb{R}^N$. La courbe de niveau est utilisée pour appliquer le seuillage et obtenir la topologie [6, 1]. La mise à jour de la topologie comprend deux étapes. La première est l'évolution de la courbe de niveau $\Psi(\boldsymbol{x}) \in \mathbb{R}^N$ par une méthode de gradient :

$$ \Psi^{(n+1)}(\boldsymbol{x}) = \Psi^{(n)}(\boldsymbol{x}) - \mathbf{P}^{(n)} \nabla \mathcal{J}^{(n)}(\boldsymbol{x}). \qquad (2) $$

Idéalement, $\mathbf{P}$ est une bonne approximation de l'inverse de la hessienne, $\mathbf{P} \approx \tilde{\mathbf{H}}^{-1}$, et $\nabla \mathcal{J}$ est le gradient (spatial) du coût. La deuxième étape est la mise à jour elle-même de la topologie. Dans le cadre de cette étude, nous employons l'application continue et différentiable :

$$ \alpha(\boldsymbol{x}) = \frac{1}{2} + \frac{1}{\pi} \arctan \frac{\Psi(\boldsymbol{x})}{\epsilon}. \qquad (3) $$

Le paramètre $\epsilon \to 0$ permet une discrimination claire des deux milieux ; ceci est utile lors du calcul de la fonction de coût à partir de la résolution du modèle physique. Au contraire, $\epsilon \neq 0$ donne une version différentiable de $\alpha(\boldsymbol{x})$, utile dans le calcul du gradient du coût $\nabla \mathcal{J}$. Dans le cadre de cette étude, nous avons donc choisi deux valeurs de $\epsilon$, l'une pour la résolution des problèmes direct et adjoint, l'autre dans le calcul du gradient. Le principal avantage est l'absence de biais important dans les modèles, et donc dans le calcul de la fonction de coût, le principale désavantage étant l'existence, d'un point de vue théorique, d'un biais entre le coût et son gradient associé : le gradient peut théoriquement ne pas être nul au minimum du coût, et le gradient peut théoriquement être nul là où le coût n'est pas exactement au minimum. Notons que l'emploi de méthodes SIMP (Solid Isotropic Material with Penalisation) aurait permis de ne pas avoir ce biais entre la fonctionnelle de coût et son gradient, mais au prix d'un biais dans la résolution du problème physique lui-même. Ceci constitue la principale raison pour laquelle une méthode basée sur les courbes de niveaux a été préférée.

Dans le cadre de cette étude, une méthode mésoscopique de lattice Boltzmann (LBM) est utilisée. La principale raison est l'application des conditions aux limites de type "bounce-back". Ainsi, contrairement aux schémas éléments finis, ou volumes finis, qui nécessitent l'emploi d'un remailleur à chaque itération du schéma itératif d'optimisation, la LBM peut, dans une certaine mesure, se passer de tels remailleurs ; le coût en temps de calcul s'en trouve abaissé. Les schémas LBM permettent aussi, assez facilement, de traiter de problèmes fluidiques multi-phasiques, plus difficiles à traiter avec des modèles macroscopiques.
Un dernier point est le calcul du gradient du coût. Pour ce faire, nous avons dérivé les équations adjointes sur les équations de Boltzmann dans leur version semi-continue : continue en espace, mais discrète selon les directions. Ceci permet, lors de l'écriture des états adjoints, d'utiliser les formes aussi compactes que possible pour tous les opérateurs (dérivée temporelle, advection, réaction) tout en utilisant facilement les opérateurs de rebond (bounce-back) là où le glissement est nul, d'entrée, et de sortie.

La physique considérée comprend un champ de température ou de concentration advecté dans un ouvert bi-dimensionnel. On se place à l'échelle mésoscopique. Pour la partie mécanique des fluides, on utilise un opérateur de collision à temps de relaxation multiples dans le schéma de lattice Boltzmann, alors qu'un schéma à temps de relaxation unique est considéré suffisant pour le champ advecté (température ou concentration). On résoudra, pour $x \in D \subset \mathbb{R}^2, t > 0$ [2] :

$$ \frac{\partial f_i}{\partial t} + c_i \cdot \nabla f_i = -M^{-1}\mathbf{S}\mathbf{M}(\mathbf{f} - \mathbf{f}^{\text{eq}}) \quad \forall i \in [1; I]. \qquad (4) $$

avec une discrétisation de l'espace des vitesses selon une quadrature du type D2Q9. Dans cette relation, $M$ est la matrice de transformation des moments, $S$ est une matrice diagonale qui comprends divers temps de relaxation, et $f^{\text{eq}}$ est la distribution à l'équilibre. Une équation similaire doit être résolue pour le champ advecté, disons en $g$. Les conditions de réflection de type rebond sont appliquées en paroi pour imposer une vitesse nulle ; d'autres conditions pour l'entrée, la sortie, pour l'état imposé ou le flux d'état imposé sont selon [7, 5]. La solution stationnaire de (4) est recherchée sur une discrétisation spatiale de type lattice, à l'aide d'un schéma discrétisé en temps. Le temps final, disons $\bar{t}$, est celui pour lequel la convergence est atteinte, c'est-à-dire lorsque $\max_x \max_i \partial_t f_i(x) \le \epsilon$.

La résolution numerique de ce modèle permet le calcul de la fonction de coût. Celle-ci, qui représente un objectif physique, adjoint à des contraintes venant du monde de l'ingéniérie, se doit d'être traduite en terme de variables d'état, en termes de $f$ et $g$ introduits plus haut, au temps final $\bar{t}$. À titre d'exemple, si l'on souhaite minimiser le flux de concentration en sortie de domaine, on cherchera $\min_{\int_{\partial D_{\text{out}}} u \cdot n c_X} d\mathbf{x}$. En variables de Boltzmann, ceci se ré-écrit en :

$$ \mathcal{J}(\mathbf{f}, g) = \frac{1}{|\partial D_{\text{out}}|} \int_{\partial D_{\text{out}}} \sum_{c_i \cdot n > 0} c_i f_i \cdot n \sum_i g_i d\mathbf{x}. \qquad (5) $$

Le gradient du coût, $\nabla \mathcal{J}$ est obtenu après résolution de la version adjointe de l'état. Comparé à (4), cet adjoint fait intervenir un signe négatif devant les opérateurs dérivée temporelle et d'advection. Ceci provient de diverses intégrations par parties [4]. Cette équation adjointe est donc résolue en temps rétrograde, à partir de $\bar{t}$, temps pour lequel la fonction de coût est calculée. Un terme source additionnel comprend les termes de couplage avec la variable de Boltzmann, $g$, qui porte le scalaire advecté, ainsi que la dérivée du coût par rapport à l'état, en l'occurence $f$ pour cette équation.

# 2 Exemples d'optimisation topologique en fluidique

Deux exemples d'optimisation topologique sont ici présentés. Le premier exemple repose sur un couplage convecto-thermique, tandis que le deuxième exemple concerne un écoulement fluide réactif.

## 2.1 Convection thermique

Dans cette première application académique, un champ de température est advecté. Plus précisément, on considère une cavité avec une entrée et une sortie. En entrée, $\partial D_{\text{in}} \subset \partial D$, la vitesse est supposée
connue, avec un profil de vitesse parabolique d'amplitude maximale donnée et telle que le nombre de Reynolds, calculé en entrée, soit égal à 20. En sortie, $\partial D_{\text{out}} \subset \partial D \setminus \partial D_{\text{in}}$, la condition est à pression nulle imposée. Ailleurs, à deux endroits distincts en frontières, $\partial D_{\text{ch}} \subset \partial D \setminus (\partial D_{\text{in}} \cup \partial D_{\text{in}})$, deux pastilles chauffantes sont maintenues à température imposée. Au contact de ces frontières, le fluide s'échauffe.

Il s'agit dans ce problème d'optimisation topologique de minimiser la température moyenne dans le domaine. Ajouté à cela, une contrainte de perte de charge est également considérée. La topologie initiale est la cavité totalement fluide. La figure 1 présente le résultat d'optimisation topologique. Au fur et à mesure des itérations, du solide est placé dans le domaine, de telle sorte que l'écoulement du fluide vienne « lécher » les parois chaudes à vitesse plus importante, et laisser en conséquence un plus grand espace sur lesquel les vitesses sont faibles et où il n'y a pas d'échauffement. Ainsi, au fur et à mesure des itérations, la fonction de coût – la température moyenne – baisse jusqu'à se stabiliser. Cela est représenté en haut de la figure, en trait plein bleu. Cette figure montre que l'optimisation topologique a permis une baisse significative du coût, augmentant en fait l'efficacité énergétique du système. D'autres exemples d'optimisation topologiques, à différents nombres de Reynolds, ainsi que des comparaisons avec des tests issus de la littérature, sont présentés dans [4].

## 2.2 Système convecto-réactif

La deuxième application concerne l'optimisation topologique de problèmes réactifs. En ce qui concerne cette étude, l'objectif futur est la conception optimale du medium de distribution des réactifs servant de support à la membrane de pile à combustible au cœur de laquelle une équation de réaction a lieu. Comme dans l'exemple précédent, on considère une cavité fluide avec une entrée et une sortie. Le nombre de Reynolds en entrée est donné et vaut 1000. La cinétique de réaction suit une loi exponentielle ; la concentration est maximale en entrée. Pour ce problème d'optimisation topologique, il s'agit de maximiser l'opération de réaction à l'intérieur du domaine. La fonction de coût associée mesure ainsi le produit de la vitesse par la concentration en sortie. On cherche donc à minimiser (5).

La figure 2 présente le résultat d'optimisation topologique. En haut à gauche sont présentés les évolutions de la fonction de coût, décroissante, et de la perte de charge. À droite, sont présentées les évolutions des topologies aux itérations 0, 50, 150, et 250. La fonction de coût a baissé d'un facteur 10 à peu près. À convergence, la topologie optimale est composée de plusieurs parties solides, localisée près de la diagonale du domaine, face au courant principal du fluide. Ainsi, le fluide s'écoule préférentiellement vers les coins diamétralement opposés, mais une partie du fluide s'écoule également au milieu du domaine, la partie solide étant perméable. Avec cette topologie optimale, la concentration en sortie de domaine est quasi-nulle. En ce sens, la réaction chimique a été drastiquement augmentée, l'efficacité énergétique du système s'en est trouvée accrue. D'autres résultats intéressant d'optimisation topologique de systèmes réactifs peuvent être trouvés dans [3].

## 3 Conclusion

Ce court papier de conférence à montré deux résultats d'optimisation topologique en fluidique. D'un point de vue formel, ces résultats reposent sur un schéma de lattice Boltzmann pour la résolution des phénomènes physiques, à l'échelle mésoscopique, sur l'introduction d'une courbe de niveaux (level-set) pour exprimer la topologie à partir d'une fonction continue et dérivable, sur la formulation adéquate du problème d'optimisation sous contraintes, et enfin sur le calcul du gradient par états adjoints. Dans les deux exemples considérés, les topologies obtenues ont permis une forte diminution de la fonction de
<figure><img src="image_3.png" /><figcaption>FIGURE 1 – Convection thermique. En haut à gauche : évolution de la fonction de coût et des pertes de charge ; à droite : évolution de la topologie aux itérations 0, 50, 100, et 250. Au milieu à gauche : module du champ de vitesse à l'itération 0 (cavité totalement fluide) ; à droite : module du champ de vitesse avec la topologie optimale. En bas à gauche : champ de température à l'itération 0 (cavité totalement fluide) ; à droite : champ de température avec la topologie optimale.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 2 – Convection-réaction. En haut à gauche : évolution de la fonction de coût et des pertes de charge ; à droite : évolution de la topologie aux itérations 0, 50, 150 et 250. Au milieu à gauche : module du champ de vitesse à l'itération 0 ; à droite : module du champ de vitesse avec la topologie optimale. En bas à gauche : champ de concentration à l'itération 0 ; à droite : champ de concentration avec la topologie optimale.</figcaption></figure>
coût, exprimant ainsi l'obtention d'une meilleure efficacité énergétique. Notons que les topologies obtenues n'auraient pu être pressenties. D'autres études nous ont permis de mettre en évidence l'existence de minima locaux. Aussi, l'implémentation d'algorithmes multi-échelles est envisagée pour convexifier la fonctionnelle à minimiser. D'un point de vue applicatif, les développements futurs concernent le multi-phasique en multi-composants, pour optimiser les canaux d'irrigation de piles à combustibles.

Références

[1] G. Allaire, F. Fe Gournay, F. Jouve, and A. Toader. Structural optimization using topological and shape sensitivity via a level set method. *Control and cybernetics*, 34(October 2004) :59, 2005.

[2] D. d'Humieres. Multiple-relaxation-time lattice boltzmann models in three dimensions. *Philosophical Transactions of the Royal Society of London. Series A : Mathematical, Physical and Engineering Sciences*, 360(1792) :437-451, 2002.

[3] F. Dugast, Y. Favennec, and C. Josset. Reactive fluid flow topology optimization with the multi-relaxation time lattice boltzmann method. submitted.

[4] F. Dugast, Y. Favennec, C. Josset, Y. Fan, and L. Luo. Topology optimization of thermal fluid flows with an adjoint lattice boltzmann method. *Journal of Computational Physics*, 365 :376 – 404, 2018.

[5] T. Inamuro, M. Yoshino, H. Inoue, R. Mizuno, and F. Ogino. A lattice Boltzmann method for a binary miscible fluid mixture and its application to a heat transfer problem. *Journal of Computational Physics*, 179 :201–215, 2002.

[6] J. A. Sethian and P. Smereka. Level set methods for fluid interfaces. *Annual review of fluid mechanics*, 35(1) :341–372, 2003.

[7] Q. Zou and X. He. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. *Physics of Fluids*, 9 :1591, 1997.