Keywords: Topology optimization, density method, bifluid, laminar flow, heat exchangers

List of figures

Chapitre 1: Introduction Générale

Motivation et Objectifs

Un grand nombre d'études et de recherches sont menées pour améliorer la performance des systèmes énergétiques et optimiser leur efficacité, en profitant de l'avancement des algorithmes mathématique d'optimisation et de la capacité de calcul des ordinateurs. Cependant 80% de l'utilisation de l'énergie est impliquée dans le transfert de chaleur. Cela souligne l'importance majeure des échangeurs de chaleur et leur implication dans presque tous les systèmes énergétiques. Ils sont utilisés dans des applications telles que les procédés, la production et la conversion d'énergie, le transport, la climatisation et la réfrigération, la récupération de chaleur et les industries de stockage et de fabrication. Les échangeurs de chaleur sont des dispositifs utilisés pour transférer de la chaleur entre deux fluides ou plus, ou entre un objet solide et un fluide, généralement sans interaction de travail.

Le plus simple échangeur de chaleur est celui à double tuyau, un pour le flux froid et un autre pour le flux chaud séparés par une certaine épaisseur. Ce type d'échangeur est le moins efficace et compact, par conséquent, plusieurs technologies de construction d'échangeurs de chaleur ont été développées au cours des années, parmi lesquelles on peut citer les plus importantes : échangeur à plaques, échangeur tubes ailettes, échangeur tubes calandre, etc.

De nombreuses études et recherches ont été menées pour améliorer le transfert de chaleur dans les échangeurs de chaleur, par l'insertion de dispositifs de perturbation de l'écoulement (ailettes persiennes, générateurs de tourbillons, etc.) et la modification de la rugosité des surfaces de transfert de chaleur. Néanmoins, il existe beaucoup de théories d'idéalisation dans la conception des échangeurs de chaleur pour obtenir une meilleure performance globale. Par exemple, concevoir des échangeurs de chaleur ayant des débits égaux dans tous les canaux [START_REF] Rohsenow | Handbook of Heat Transfer[END_REF], ou régler l'équilibre du coefficient de transfert thermique de chaque côté de l'échangeur de chaleur [START_REF] Rao | Effect of flow distribution to the channels on the thermal performance of a plate heat exchanger[END_REF]. Cependant, l'amélioration de la performance des échangeurs de chaleur est une procédure plus sophistiquée que de simplement utiliser des techniques pour augmenter le coefficient de transfert de chaleur, ou considérer des théories d'idéalisation. Ceci est dû à de nombreux facteurs intervenant dans l'échangeur de chaleur. Le taux de transfert de chaleur et la perte de charge due au frottement du fluide, sont les phénomènes physiques, toujours en opposition, influençant le plus les caractéristiques de l'échangeur. Ainsi l'optimisation des échangeurs de chaleur a été un domaine d'étude et de recherche intensif, pour sa capacité à améliorer la performance des échangeurs de chaleur et leur efficacité en prenant compte de tous les facteurs en opposition et les limitations de conception. C'est pourquoi l'amélioration des performances des systèmes, y compris les échangeurs de chaleur, dépend de la capacité à répondre aux spécifications demandées pour l'échangeur de chaleur lui-même. L'objectif de cette recherche est de repousser les limites des connaissances et de développer des outils et des méthodes pour permettre la création d'une nouvelle génération d'échangeurs de chaleur. Le nouveau concept développé est basé sur les méthodes les plus récentes et les plus complexes dans le domaine de l'optimisation de la configuration, les techniques d'optimisation topologiques qui ne sont pas basées sur une géométrie prédéfinie. Ces méthodes permettent d'atteindre une architecture complexe et efficace basée strictement sur les objectifs et contraintes définis. L'état actuel des travaux scientifiques permet l'application de l'optimisation topologique aux échangeurs de chaleur comprenant un seul fluide et un solide. Le présent travail vise à étendre les méthodes d'optimisation topologique en mécanique des fluides à des cas incluant deux fluides, ce qui est le cas pour les échangeurs de chaleur fluide-fluide. Cependant, une attention devrait être accordée à la complexité des géométries générées par l'optimisation de la topologie. Une grande question se pose donc: comment ces structures seront-elles fabriquées en particulier pour la production à grande échelle? L'avancement dans la technologie de fabrication additive est la réponse évidente à cette question. Il est donc très utile d'associer l'impression 3D à l'optimisation topologique, pour le développement d'une nouvelle méthode innovante dans la conception et l'optimisation des échangeurs de chaleur.

Optimisation des échangeurs de chaleur

L'optimisation est le mécanisme de sélection de la meilleure solution dans une situation particulière soumise à un certain nombre d'obstacles et de limitations. Le critère qui définit la meilleure solution est la fonction objectif. Les limitations sur les solutions disponibles sont définies par les contraintes. Ce qui décrit différentes solutions sont les variables de problème auxquelles nous essayons d'assigner les meilleures valeurs pour minimiser ou maximiser la fonction objectif. La façon dont nous pouvons réaliser le processus d'optimisation est définie par l'algorithme d'optimisation que nous utilisons. Si l'on veut optimiser une fonction objectif f (x) = x, la meilleure solution est simplement infinie. Ainsi, un problème d'optimisation n'a aucun sens s'il n'y a pas de conflit entre plusieurs fonctions objectifs ou entre une fonction objectif et une contrainte. De même, dans les échangeurs de chaleur, s'il n'y a pas de limitation sur le volume ou la masse ou si l'on ne tient pas compte de la perte de charge, l'échangeur de chaleur optimal pour avoir un transfert de chaleur maximal est celui ayant une longueur infinie. Les variables du problème d'optimisation des échangeurs de chaleur peuvent être les conditions de fonctionnement de l'échangeur, propriétés physiques des matériaux et les fluides et les paramètres géométriques du dispositif.

Les problèmes d'optimisation, dans lesquels les paramètres géométriques sont les paramètres d'optimisation, sont classés en trois catégories selon le degré de liberté et la possibilité de modifier la géométrie: optimisation de taille, de forme et topologique. Dans l'optimisation de la taille, les variables du problème mathématique sont les paramètres géométriques de la structure telle que la longueur, la largeur, le rayon, etc. La forme et la connectivité des éléments entre elles sont connues et fixées. Par conséquent, la solution optimale finale est très similaire à la conception de base initiale. L'optimisation de forme augmente le degré de liberté du problème, où elle peut changer la taille et la forme simultanément en ajoutant des variables capables de déformer la forme de la structure (par exemple modifier la forme des canaux d'écoulement). Cependant, l'architecture de la structure est toujours similaire à la conception initiale, puisque la topologie globale de l'échangeur de chaleur ne peut pas être modifiée. Enfin, en optimisation topologique, chaque maille du domaine d'optimisation est un paramètre de conception; ce qui permet d'ajouter ou de supprimer du matériel dans chaque point de l'espace de conception sans être limité à une topologie initiale. Cela augmente considérablement le nombre de variables dans le problème d'optimisation, ce qui rend la convergence plus difficile. Cependant, l'optimisation topologique est devenue si attrayante pour sa capacité à atteindre des configurations innovantes et complexes basées strictement sur les objectifs et contraintes définis. L'optimisation de la taille et de la forme a été largement appliquée à l'optimisation des échangeurs de chaleur, tandis que l'optimisation topologique est encore limitée aux échangeurs de chaleur fluide-solide.

Les deux critères principaux les plus utilisés dans les problèmes d'optimisation des échangeurs de chaleur sont les suivants : les critères basés sur les irréversibilités thermodynamiques, comme la minimisation de la génération d'entropie et/ou de la dissipation de l'entrasse, et les critères économiques qui consistent à minimiser le coût du fonctionnement et le coût initial de l'échangeur. Il existe de nombreux autres critères d'optimisation, comme la diminution de la perte de charge totale, augmentation du taux de transfert thermique, diminution de la température maximale ou moyenne, ou par exemple minimisation de la masse totale de l'échangeur de chaleur dans l'industrie aéronautique, etc.

Optimisation topologique

L'optimisation topologique a été développée à l'origine pour l'optimisation des problèmes de structures mécaniques. L'objectif était de trouver la forme qui utilise le minimum de matière tout en maintenant les contraintes mécaniques inférieures à un niveau acceptable. L'optimisation topologique a été définie par Bendsoe et Sigmund [START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF] comme une optimisation de forme des structures qu'elle devrait définir à chaque point de l'espace de conception s'il existe un matériau ou non, la topologie de la structure n'étant pas fixée a priori. Ainsi, à partir d'un domaine initial vide, complet ou dans n'importe quel état intermédiaire, les paramètres de contrôle utilisés permettent de créer sans restriction des créations de trous et d'agglomérats de matériaux, afin de trouver la meilleure topologie possible. Récemment, le concept d'optimisation topologique a été appliqué à un large éventail de disciplines physiques comme les mécaniques des fluides, le transfert de chaleur, l'acoustique, l'électromagnétique et l'optique. Dans les problèmes de transport de masse et chaleur en optimisation topologique, l'objectif était de trouver l'architecture optimale qui correspond à un compromis entre la minimisation de la perte de charge et la maximisation de transfert de chaleur. Cependant, la mise en oeuvre de méthodes d'optimisation topologique est assez complexe car elle nécessite un algorithme capable d'allouer et de réallouer efficacement le matériel dans un domaine ayant les dimensions et les conditions aux limites prédéfinies.

En optimisation topologique, l'espace de domaine est discrétisé en des petites mailles, également appelés cellules, où chaque cellule contient une variable de conception adimensionelle. Les valeurs de toutes les variables de conception dans toutes les cellules définissent la forme de la structure entière. Ainsi, le problème d'optimisation est de trouver les valeurs optimales de toutes les variables de conception, en minimisant une certaine fonction objectif en respectant les fonctions contraintes, qui sont généralement des fonctions de porosités qui limitent le volume maximal de l'une des matières dans le domaine d'optimisation.

Il existe deux grandes familles de méthodes de la résolution du problème d'optimisation topologique: les approches discrètes et les approches continues, qui correspondent respectivement aux méthodes dans lesquelles le paramètre local contrôlant le matériau dans chaque cellule peut prendre des valeurs discrètes, ou des valeurs bornées continues. Ces méthodes dépendent du gradient des fonctions objectifs et des contraintes. Le grand nombre de variables de conception et la nécessité de calculer la dérivée totale de la fonction objectif par rapport aux variables nécessitent des techniques mathématique avancées pour assurer la convergence du vecteur de variables vers la solution optimale. L'optimisation topologique est devenue un domaine de recherche bien développé avec de nombreuses techniques pour traiter les problèmes d'instabilité numérique fréquemment rencontrés dans l'optimisation topologique tels que les damiers, les dépendances de maillage et les minima locaux.

Méthodes d'optimisation topologiques

Les méthodes d'optimisation topologiques visent à forcer le paramètre de conception à prendre progressivement des valeurs discrètes, éliminant ainsi les régions grises et conduisant à un domaine noir et blanc, où les variables de conception dans chaque cellule sont égales à 0 ou 1 (solution 0-1). Cependant, il existe aussi des méthodes qui peuvent résoudre des problèmes combinatoires discrets et qui sont appelés approches discrètes. Nous allons considérer un problème d'optimisation topologique où nous cherchons la distribution optimale des deux phases A et B. Dans les approches discrètes, la phase à l'intérieur de la cellule est modifiée en une seule étape entre A et B (hard-kill). Dans les approches continues, une petite quantité de phase A est remplacée par la phase B ou vice-versa dans chaque étape, jusqu'à ce que nous atteignions à la fin du processus d'optimisation une cellule complètement correspondante à A ou B (soft-kill).

Les méthodes d'optimisation topologique reposent sur trois parties principales: le solveur direct du problème physique (éléments finis, volumes finis, ...), la méthode d'analyse de sensibilité (adjoint discret, adjoint continu, ...) pour calculer la dérivée totale des fonctions objectifs et contraintes par rapport aux variables de conception, et un optimiseur numérique.

La méthode la plus rencontrée en optimisation topologique en transfert de masse et chaleur est la méthode de densité. Cette méthode consiste à utiliser une fonction d'interpolation pénalisée pour calculer les quantités physiques dans chaque cellule, par ex. la rigidité du matériau dans les problèmes de structure mécanique et la conductivité thermique dans les problèmes de conduction de chaleur etc., en fonction des variables de conception continue. Le principal défi des méthodes de densité est l'introduction d'une fonction d'interpolation capable d'orienter la solution vers des valeurs 0-1 discrètes et d'omettre des valeurs intermédiaires de la variable de conception, tout en assurant une représentation physique réelle des matériaux fictifs correspondant à des densités intermédiaires, connus sous le nom de matériaux gris, qui doivent être entièrement éliminés quand le problème converge vers la solution finale. Un schéma d'interpolation populaire pour satisfaire les conditions ci-dessus est la formule SIMP (Solid Isotropic Material with penalization).

Les méthodes de level set sont des techniques computationnelles introduites en 1988 par Osher et Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] pour le suivi des interfaces mobiles. L'idée principale des méthodes level set est d'introduire une fonction dépendant du temps et de l'espace qui définit l'interface entre les deux matériaux présents dans le problème d'optimisation.

L'approche évolutive ESO (Evolutionary structural optimization) qui utilise des variables de conception discrètes, a été introduite par Xie et Steven [START_REF] Xie | A simple evolutionary procedure for structural optimization[END_REF] pour l'optimisation des structures mécaniques. Cette méthode est basée sur le concept simple de retrait progressif d'un matériau inefficace d'une structure jusqu'à ce que la contrainte déterminant le volume de matériau dans le domaine de conception soit satisfaite. Yang et al. a développé l'Optimisation structurelle évolutive bidirectionnelle (BESO), une version améliorée de l'ESO, qui permet de retirer et d'ajouter le matériau simultanément [START_REF] Young | 3D and multiple load case bidirectional evolutionary structural optimization (BESO)[END_REF]. Le retrait et l'addition de matière sont basés sur la valeur du nombre de sensibilité. Les approches évolutives (ESP et BESO) ont été largement appliquées aux problèmes de structures mécaniques avec de nombreuses techniques développées pour faire face aux difficultés rencontrées en raison de l'aspect discret du problème. Dans la conduction de chaleur pure, les approches évolutives étaient également applicables mais les résultats ont montré que la méthode conduit à l'optimum local. Par conséquent, les approches évolutives ne se sont pas intéressantes pour l'optimisation de la topologie de la mécanique des fluides. Malgré l'attention portée aux techniques évolutives dans les problèmes de structure mécanique et problèmes de transfert de chaleur par conduction, elles n'ont pas été prises en compte dans les problèmes d'écoulement des fluides selon la revue de la littérature. La méthode Level Set a été trouvée attrayante pour les problèmes d'écoulement en raison de leurs résultats dans des simulations numériques 2D et 3D pour divers types de flux. Cependant, cette méthode peut seulement évoluer à partir des interfaces existantes et n'est pas capable de générer de nouveaux trous, ce qui signifie qu'il est impossible de générer de nouveaux canaux dans l'optimisation des écoulements. Ceci est considéré comme un inconvénient conceptuel de la méthode, surtout si elle sera utilisée pour l'optimisation topologique des échangeurs de chaleur. La nucléation de nouveaux trous dans Level Set a été possible en la combinant avec la méthode Topological Sensitivity. Cette méthode combinée a été appliquée et testée par divers auteurs. Les résultats montrent que la solution finale reste fortement dépendante de l'estimation initiale.

Conclusion

La méthode de densité est complètement indépendante de l'estimation initiale et la génération de canaux et de structures complexes dépend uniquement des fonctions d'objectifs et des contraintes. De plus, la revue littérature a montré que la méthode de densité a été appliquée sur la majorité des problèmes liés au transfert de chaleur et de masse. Malgré la nécessité d'un temps de calcul élevé, la méthode de densité semble être la méthode la plus appropriée pour étendre l'application de l'optimisation topologique en mécanique des fluides au domaine bifluide, ce qui n'était pas envisagé auparavant.

Plan de la thèse

L'algorithme général de la méthode de densité est composé de trois étapes principales :

 Le solveur CFD utilisant la méthode des volumes finis.  L'analyse de sensibilité basée sur la méthode d'adjoint discret.  La méthode des asymptotes mobiles comme optimiseur numérique.

Le reste du document est divisé comme suit:  L'algorithme détaillé de la méthode d'optimisation présenté ci-dessus et le développement détaillé de chaque partie de la méthode sont présentés au chapitre 2. Deux formulations différentes seront comparées, l'une utilisant une variable de conception unique dans chaque cellule de conception et la seconde utilisant deux variables de conception. Dans chaque cellule de conception qui double le nombre de variables du problème.  La séparation des fluides sera examinée au chapitre 3, où chaque fluide doit prendre son propre trajet dans le domaine d'optimisation indépendamment de l'autre.  Dans le chapitre 4, la maximisation du transfert de chaleur entre les deux fluides séparés sera considérée.  Enfin conclusion et perspectives dans le chapitre 5.

Chapter 1

General Introduction

1.1. Introduction

Motivation

In the 2015 United Nations climate change conference held in Paris, representatives of 196 nations adopted a long term strategy to respond to the threats of climate change and deal with greenhouse gas emissions mitigation plans. The agreement, known as "Paris agreement", aimed to limit the rise in global average temperature by holding it this century below 2°C above preindustrials level [START_REF]The Paris Agreement -main page[END_REF]. As part of this agreement, the French environment minister announced in July 2017 his country plan to attain neutral carbon equilibrium in 2050, by reducing human carbon emissions to the level of ecosystems carbon's absorption capacity [START_REF]Plan Climat[END_REF]. The French plan also considered a four billion Euros investment to increase energy efficiency and stop coal usage for electricity production by 2022 [START_REF]The Independent[END_REF]. Beside the problems related to global warming and climate change, energy management policies face various challenges. First, population growth increases the demand for energy services [4]. Furthermore, the increase in the ratio of urban population to rural population augments the demand on energy even more. Second, an increase in gross domestic product (GDP) is associated with an increase in energy consumption, which tend to vary according to the GDP growth in different economy sectors [START_REF]orld Energy Outlook 2016[END_REF]. Energy market is also influenced by many other sectors, as technology innovations, oil and gas prices, carbon emissions pricing by some governments, etc. All these reasons explain the studies and researches conducted to improve the performance of energy systems and optimize their efficiency, by taking advantage of advancements in mathematical optimization tools and computers calculation capacities.

On the other hand, about 80% of energy utilization is involved in heat transfer. This highlights the major importance of heat exchangers and their involvement in nearly every energy system. They are used in applications such as processes, energy production and conversion, transport, air conditioning and refrigeration, heat recuperation and storage and manufacturing industries. Heat exchangers are devices used to transfer heat between two or more fluids, or between a solid object and a fluid usually without work interactions. Usually, fluids do not mix in heat exchangers, and heat is transferred through a dividing wall without fluid leakage. However, there are still some types of heat exchangers where fluids enter in direct contact and are later separated. In this case, heat transfer is mainly caused by phase change enthalpy [START_REF] Shah | Fundamental of heat exchanger deisng[END_REF]. Figure 1.1 shows a general representation of a fluid to fluid heat exchanger. and stand respectively for temperature, pressure and thermal heat transfer load. Subscripts and stands respectively for inlet, outlet, cold and hot. The heat exchanger is characterized by the total heat thermal power transferred from the hot stream to the cold stream, the pressure drop of both fluids: + , and other geometrical parameters like its mass, volume and compactness. The simplest heat exchanger is a double pipe exchanger, with one pipe for cold stream and another one for hot stream. This type of heat exchangers has the lowest efficiency and compactness, hence, a wide range of heat exchangers construction technologies have been developed over the years. Among these technologies, the following are the most known and used (Figure 1.2): tubular heat exchangers, as shell and tubes used for high pressure and temperature flow conditions [START_REF] Shah | Fundamental of heat exchanger deisng[END_REF], plate type heat exchangers characterized by a high transfer coefficient but cannot endure high pressure and temperatures flows neither high temperature gradient. There exist many other technologies as extended surface heat exchangers that use fins to increase heat transfer surface, regenerators and adiabatic wheels where heat transfer process is not continuous, etc. Heat exchangers are also characterized by their flow arrangement. We can distinguish three different types: parallel flow where the fluids flow parallel to each other in the same direction, counter flow where fluids flow parallel but in opposite direction and cross flow where fluids flow in perpendicular direction. Heat exchanger design is a complex iterative procedure, due to many physical phenomena occurring inside the device, and due to the interactions between these phenomena and their mutual dependency. Thermal design of heat exchangers aims to determine the required heat transfer surface for a fixed heat load duty, or determine the rate of heat transfer for a fixed heat transfer surface. The most used basic thermal design methods are logarithmic mean temperature difference method (LMTD) and effectiveness-number of transfer units' method (ε-NTU). Mechanical design has also a major importance in heat exchangers design. It aims to handle thermal and pressure stresses, and ensure durability of the device at different operational phases [START_REF] Shah | Fundamental of heat exchanger deisng[END_REF].

Many studies and researches were conducted to enhance heat transfer in heat exchangers, by insertion of flow disturbance devices and modification of the roughness of heat transfer surfaces. These techniques enhance heat transfer by making the flow turbulent near heat transfer surface, by breaking the laminar layer of the flow to reduce the thermal resistance and by increasing the residence time of heat transfer fluids. Among these devices we mention: vortex generators [START_REF] Huisseune | Performance enhancement of a louvered fin heat exchanger by using delta winglet vortex generators[END_REF] (Figure 1.3), louvered fins [START_REF] Carija | Heat transfer analysis of fin-and-tube heat exchangers with flat and louvered fin geometries[END_REF] twisted tapes [START_REF] Liu | A comprehensive review on passive heat transfer enhancements in pipe exchangers[END_REF], ribs [START_REF] Gawande | Effect of roughness geometries on heat transfer enhancement in solar thermal systems -A review[END_REF] [START_REF] Taslim | Experimental Heat Transfer and Friction in Channels Roughened With Angled, V-Shaped and Discrete Ribs on Two Opposite Walls[END_REF], spiral fins [START_REF] Pongsoi | Heat transfer and flow characteristics of spiral fin-and-tube heat exchangers: A review[END_REF], circular fins [START_REF] Jang | The thermal-hydraulic characteristics of staggered circular finned-tube heat exchangers under dry and dehumidifying conditions[END_REF], etc. Nevertheless, there exist a lot of idealization theories in heat exchangers design to get a better overall performance. For example designing heat exchangers with equal flow rates in all channels [START_REF] Rohsenow | Handbook of Heat Transfer[END_REF], or setting equilibrium in heat transfer coefficient at each side of the heat exchanger [START_REF] Rao | Effect of flow distribution to the channels on the thermal performance of a plate heat exchanger[END_REF], etc. Enhancing the performance of heat exchangers is quite a more sophisticated procedure than simply using techniques for increasing heat transfer coefficient, or considering idealization theories. This is caused by many trade-off factors occurring inside the heat exchanger. The most important one is the trade-off between heat transfer rate and pressure drop due to fluid friction. If the designer wants to increase heat transfer rate, he could simply decrease the diameter of the tubes but the pressure drop will then increase. In this case, what determines the enhancement of the heat exchanger's performance is the cost of the pumping power due to the pressure drop versus the benefit of the recovered thermal power. Another example is the trade-off between the capital cost and the operating cost where a more efficient heat exchanger may decreases the operating cost but requires a higher capital cost and vice versa. Design of heat exchangers is also subject to many constraints, which vary according to the application, like total weight and volume, design limitations to avoid corrosion, fatigue failure, etc. Optimization of heat exchangers consists of finding a compromise between all trade-off factors within the feasible solutions that respect all design constraints, by minimizing a certain objective or optimization criteria.

Heat exchangers optimization has been an intensive field of study and research, for its capability to improve the performance of heat exchangers and their efficiency while taking into account all trade off factors and design limitations. Dimitrios et al. [START_REF] Misirlis | Optimization of Heat Exchangers for Intercooled Recuperated Aero Engines[END_REF] conducted an optimization of heat exchangers mounted on the hot gas exhaust nozzle of an intercooled recuperated aero engine. The optimization resulted in two new recuperators, which were compared with the initial baseline design on the basis of their weight and specific fuel consumption of the aero engine. The initial non-optimized heat recuperator was capable of achieving 12.3% reduction in fuel consumption in relation to a non intercooled aero engine. The first optimized recuperator increased fuel consumption reduction to 13.1% in relation to non intercooled aero engine, while reducing the weight of the recuperator by 5%. The second optimized recuperator was less efficient regarding fuel consumption, whose reduction in relation to a non intercooled engine dropped to 9.1%, but on the other hand the total weight was reduced by 50% in relation to the initial non optimized recuperator. Ghadamian et al [18] optimized the operation conditions of heat exchanger used for heat recuperation in a cement industry. They were able to increase heat recuperation by 592.2 Kw/year without any increase in cost. Heat exchanger used for waste heat recovery in industry was also considered by Yildirim and Soylemez [19]. Their resulted optimized plate heat exchanger achieved approximately 0.91 million $ net profit over 10 years life cycle, whereas the initial heat exchanger was only capable of saving around 0.78 millions.

Caputo et al. [START_REF] Caputo | Heat Exchanger Optimized Design Compared With Installed Industrial Solutions[END_REF] tested the effectiveness of an optimization method on several heat exchangers installed in chemical plants. A shell and tube Soda-water heat exchanger weight was reduced from 5.287 Kg to 2.697 Kg, whereas pressure drop was decreased from 5.5 kPa to 1.7 kPa on shell side and remained approximately the same on tube side. Another potassium hydroxide-water heat exchanger was studied, where a 58% reduction in weight was achieved. The pressure drop in the optimized heat exchanger was increased from 5 to 8 kPa on the shell side, but on the tube side the pressure drop was significantly reduced from 48.4 kPa to 7 kPa. All this improvement resulted in a better performance regarding operational cost, and the heat exchanger became shorter and more compact. [START_REF]Design and multi-objective optimization of heat exchangers for refrigerators -ScienceDirect[END_REF] provided a multi-objective optimization of a heat exchanger used in refrigeration where they showed the trade-off between energy consumption and material cost. Regarding energy consumption, the best design presented a reduction of 8.92 % in relation to the baseline design, but it needed a 50.19 % increase in material cost. The best achievement in term of material cost was a reduction by 41.82% at the expense of a 6.15% increase in energy consumption on a daily basis. In that case, heat exchangers optimization provides best trade-off solutions, and the choice of a final design is based on the designer strategy regarding the competing objectives. This brief literature review shows the advantage and profit gained by using advanced optimization tools in heat exchangers design.

Gholap and Khan

Research objectives

As seen in the last paragraph, improving the performances of systems including heat exchangers depends on the ability to meet the specifications requested for the heat exchanger itself. The objective of this research is to push the limits of the knowledge base and develop tools and methods to enable the creation of a new generation of heat exchangers. The new design concept developed is based on most recent and complex methods in the field of configuration optimization, the topology optimization techniques which are not based on predefined geometry. These methods allow reaching a complex and efficient design based strictly on the defined objectives and constraints. The current state of scientific work allows the application of topology optimization to heat exchangers including a single fluid and a solid. The present work aims to extend the topology optimization methods in fluid mechanics to cases including two fluids, which is the case for fluid-fluid heat exchangers. Meanwhile, an attention should be given to the complexity of geometries generated by topology optimization. Hence, a big question arises: how these structures will be manufactured especially when it comes to large scale production? The advancement in additive manufacturing technology is the obvious answer to this question. Additive manufacturing is the process of adding layer upon layer of a given material. It reads information from a computer-aided design (CAD) file to add successive layers of materials to fabricate the designed object. The first use of additive manufacturing was limited to create prototype or visualize a part for presentations purposes. However, currently the additive manufacturing is used to produce end-use products for a wide range of applications such as aircraft parts, automobile, medical equipments etc. Therefore, the advance in 3D printers technologies raised interest in topology optimization development and application in various engineering industries. It is then very useful to associate 3D printing with topology optimization, for the development of a new generation in heat exchangers design and optimization. Optimization is the mechanism of selecting the best solution in a particular situation subject to a number of obstacles and limitations. The criterion that defines the best solution is the objective function. The limitations on available solutions are defined by the constraints. What describe different solutions are the problem variables to which we are trying to assign the best values to minimize or maximize the objective function. How we can achieve the optimization process is defined by the optimization algorithm we're using. If one wants to optimize an objective function , the best solution is simply infinite. Hence, an optimization problem has no sense if there is not a conflict between many objective functions or between an objective function and a constraint. Similarly, in heat exchangers, if there is no limitation on volume and mass or there is no consideration to pressure drop, the optimal heat exchanger regarding heat transfer is the one having an infinite length. Figure 1.4 represents a general schematic for heat exchangers optimization. In next paragraphs most encountered optimization variables, criteria and numerical algorithms in heat exchangers optimization are presented.

Heat exchangers optimization 1.2.1. Introduction to optimization problem

Optimization variables

The design variables in heat exchangers optimization problems are classified as follows:  Operating conditions of the heat exchanger (mass flow rate of hot or cold stream, terminal temperatures, etc…)  Material properties (thermal conductivity, surface roughness, etc….)  Geometrical parameters that defines the optimal architecture of the heat exchanger. Optimization problems, in which the geometrical parameters are the optimization parameters, are classified under three categories according to the degree of freedom, and the capability of changing the geometry: size, shape and topology optimization, represented in Figure 1.6. In size optimization the problem variables are the geometrical parameters of the structure such as length, width, radius etc. The shape and the connectivity of the elements between them are known and fixed. Hence, the final optimal solution is very similar to the initial baseline design.

Shape optimization increases the degree of freedom of the problem, where it can change the size and the shape simultaneously by adding variables able to deform the boundaries of the structure. However, the architecture of the structure is still similar to the initial design, since the global topology of the heat exchanger cannot be varied. Finally, in topology optimization, every mesh in the optimization domain is a design parameter; which allows adding or removing material in every point in the design space without being limited to an initial topology. This increases significantly the number of variables in the optimization problem what make it more difficult to converge, as seen in Figure 1.5. However, topology optimization has become so appealing for its capacity to reach innovative and complex configurations based strictly on the defined objectives and constraints. Size and shape optimization have been widely applied on optimization of heat exchangers, whereas topology optimization is still limited to fluid to solid heat exchangers.

Numerical optimization algorithms

Optimization techniques were subject to a considerable progress since early fifties of the twelve century, with the advancement in digital computers and electronic calculation capacities. Optimization algorithms vary according to the information and the method used to search for the next better solution. There exist in literature a wide variety of available optimization algorithms, classified under many criteria, as their capability to reach local of global optimum, their dependency or no on gradient information, the use of single trajectory or a population, etc. In the following, the most used optimization algorithms in engineering problems will be presented and classified as gradient-based and gradient-free algorithms.

 Gradient-based algorithms: this family of optimization algorithm uses the derivative of the objective function as a direction to reach the optimum. The major advantage of these methods is that they can solve optimization problems with an extremely large number of variables. However, their major drawbacks are that they can locate a local optimum and that they are complex to implement. The essential part in calculation time in these algorithms is dedicated to the evaluation of the sensitivity of the objective function. Among most popular gradient based algorithms used, Fletcher-Reeves for unconstrained problem and the Sequential Linear Programming (SLP) and Sequential Quadratic Programming (SQP) [START_REF][END_REF] for constrained optimization problems..  Gradient-free algorithms: the most popular family of this type of optimization techniques are the evolutionary algorithms, which are based on phenomena from nature to evolve toward optimal solution. Evolutionary algorithms are easy to implement, and can reach usually global or near global optimal solution. On the other side, evolutionary algorithms suffer from not being able to handle high number of variables and are computational costly. Genetic algorithm [START_REF]Wiley: Evolutionary Optimization Algorithms -Dan Simon[END_REF] is the most efficient and popular type of evolutionary algorithms. It is inspired from science of genetics for the survival of the fittest, more specifically Darwin's theory and Mendel's law for genetic evolution and inheritance. It uses biological operators such as crossover, mutation and selection. Design points having the best performance regarding the objective function, are used for the next generation of the optimization. Other popular evolutionary algorithms frequently used are the Particle Swarm Optimization and simulated annealing.

Optimization criterions

In this paragraph, the objective functions usually encountered in the optimization of heat exchangers are presented. The two main criterions are thermodynamic and economic aspects, which could also be coupled in the same problem as it will be seen later. There exist many other optimization criterions depending on the application of the heat exchanger, like decreasing the total pressure drop, increasing the heat transfer rate, decreasing the maximum temperature or average temperature, or for example the minimization of total heat exchanger weight in aeronautical industry, etc.

Thermodynamic criteria

We will start by presenting the thermodynamic irreversibility occurring in a heat exchanger. In heat exchangers, there are two types of losses: losses associated to irreversibility due to heat transfer between two fluids having a certain temperature difference, and losses associated to mass transfer due to friction between fluids and internal walls of the heat exchanger. Regarding fluid flow, the mechanical energy is not conserved, and a part of it is transformed into thermal energy. This part of energy lost during fluid transport process is considered as a measure for irreversibility. Usually, in most physical phenomena, the quantity of energy turned into heat is considered as the quantity of irreversibilities. When it comes to heat transfer itself, thermodynamic quantities as entropy production and exergy destruction are considered as irreversibilities. Heat exchangers optimization using thermodynamic criteria are based on the minimization of these quantity of irreversibilities, more specifically entropy generation and entransy dissipation.

Entropy production was developed by Bejan [START_REF]Entropy generation minimization: The new thermodynamics of finite-size devices and finite-time processes[END_REF] who defined it as "the thermodynamics imperfection to heat transfer, mass transfer and fluid flow irreversibilities". Entropy production parameter takes into account both transport processes, heat transfer and mass transfer. Hence, in any variation in the heat exchanger's geometry to increase heat transfer, associated mechanical energy dissipation is simultaneously evaluated and taken into consideration. It should be noted that the distribution of entropy production itself inside the exchange influences the overall exchanger effectiveness. In fact, the minimal total entropy production of the entire system corresponds to the one having the most uniformly possible local entropy production distribution [START_REF]Equipartition of entropy production. An optimality criterion for transfer and separation processes[END_REF].

Entransy dissipation is another physical quantity developed in 2007 by Guo et al. [START_REF] Guo | Entransy-A physical quantity describing heat transfer ability[END_REF] to measure the irreversibilities in a heat transfer process. It is developed to measure the ability of an object (fluid in case of a heat exchanger application) to transfer heat in analogy with the electric capacitance of a body, which describes its charge transfer ability. Xu et al. [START_REF] Xu | Application of entransy dissipation theory in heat convection[END_REF] applied the entransy dissipation theory of Guo. et al. on internal and external fluid flow. Hence, entransy dissipation was able to take into account irreversibilites due to heat and mass transfer simultaneously.

Before the introduction of entransy dissipation, entropy generation number was considered as the main criterion in heat exchanger optimization based on thermodynamic irreversibilities minimization. The recent studies have shown a preference of methods based on entransy dissipation over entropy generation methods. Twenty different heat exchangers were analyzed by Qian et Li [START_REF] Qian | Analysis of entransy dissipation in heat exchangers[END_REF]. The results showed that the minimum entransy rate corresponds in all the heat exchangers tested to the highest heat transfer rate. The minimum entropy generation suffered for many cases from the entropy generation paradox, where the efficiency of the heat exchangers can be at its maximum, minimum or anything between when entropy generation reaches its minimum [START_REF] Shah | Entropy Generation Extrema and Their Relationship With Heat Exchanger Effectiveness-Number of Transfer Unit Behavior for Complex Flow Arrangements[END_REF]. However, it was demonstrated that the two physical quantities are needed to evaluate irreversibility in heat transfer [START_REF] Guo | Effectiveness-thermal resistance method for heat exchanger design and analysis[END_REF]. When the purpose of heat transfer is for heat-work conversion, the entropy generation is a better irreversibility measurement, whereas the entransy dissipation is better when the heat transfer is for heating and cooling purposes [START_REF] Guo | Effectiveness-thermal resistance method for heat exchanger design and analysis[END_REF] [31].

Economic criterion

The economic criterion consists of minimizing the total cost of a heat exchanger which is the sum of the capital cost and the operating cost. Economic objective is a widely used criterion in optimization of heat exchangers as it reflects directly the purpose of using heat exchangers in the majority of energy systems. In some applications, it is necessary to take into consideration the economic benefit from each degree of temperature recovered by the heat exchanger, versus the price of each kW needed to operate the exchanger. This is not possible in the methods based on the minimization of irreversibility rate alone. The cost function could be used alone as single objective function [START_REF] Swamee | Optimum design of double pipe heat exchanger[END_REF] or used in multi-objective design optimization to find a trade-off between the exchanger cost and its effectiveness [START_REF] Ghanei | Thermaleconomic multi-objective optimization of shell and tube heat exchanger using particle swarm optimization (PSO)[END_REF].

Literature review on fluid-to-fluid heat exchangers optimization

Size optimization

Huang et al. [START_REF] Huang | Optimal design of vertical ground heat exchangers by using entropy generation minimization method and genetic algorithms[END_REF] optimized a vertical ground heat exchanger used in HVAC systems, by minimizing entropy generation number using a genetic algorithm. The variables are a set of geometrical parameters and the material conductivity. The entropy generation expression is coupled with the irreversibility caused by fluid friction and heat transfer as a single objective function. His optimal design has an entropy generation number of 12.2% less than the initial design. He analyzed the advantage made by the optimization method from an economical point of view over a 10 years operation period. The results showed that the capital cost is 1.67 % higher for the optimal design, but the operation cost was 7.2% lower. Thus, he achieved a 5.5% net profit in total cost over the operation period. Guo et Xu [START_REF] Guo | The application of entransy dissipation theory in optimization design of heat exchanger[END_REF] applied theory of entransy dissipation on size optimization of a shell and tube heat exchanger using a genetic algorithm. He also showed the benefit of splitting entransy dissipation due to heat transfer and fluid flow as two objective functions and used them in a multi-objective optimization instead of a single objective optimization. The advantage of a multi-objective function is that the designer can control the preferences of maximization of heat transfer and minimization of pressure drop. Results showed that in the design of a heat exchanger with fixed heat load, the single objective optimization improves the performance of the heat exchanger. However, when the heat transfer area is fixed, the improvement of the heat exchanger effectiveness is at the expense of increasing the pumping power. The multi objective optimization design can achieve the same effectiveness as single objective design with less consumption in pumping power, in case of fixed heat transfer area.

Huang [START_REF] Huang | A multi-objective design optimization strategy for vertical ground heat exchangers[END_REF] compared two different optimization methods: a single objective optimization in which entropy generation is the optimization criterion, and a multi objective method where the entropy generation and the total heat exchanger cost are the objective functions. Optimization procedures were applied on a vertical ground heat exchanger to find the optimal values of various geometrical parameters using a genetic algorithm. The heat exchanger optimized by a single objective method has an operating cost 0.8 % lower than a heat exchanger optimized by the multi-objective method, and a 0.82% lower entropy generation number. On the other side, the capital cost is 10% lower for the heat exchanger optimized by the multi-objective method.

Hence, the multi-objective optimization design has a 6.22% lower total cost in relation to the design of the single objective optimization method, and a 9.5% lower total cost in relation to the original design. This result shows the importance of associating the economical criterion to the thermodynamic quantities in size optimization of heat exchangers. Thermal-economic multiobjective optimization was also considered by Ghanei et al. [START_REF] Ghanei | Thermaleconomic multi-objective optimization of shell and tube heat exchanger using particle swarm optimization (PSO)[END_REF] in size optimization of a shell and tube heat exchanger using particle swarm optimization algorithm. He achieved a 7.36% improvement in the total heat transfer rate while reducing the total cost by 0.36%. Juan and Qin [START_REF] Juan | Multi-objective optimization of a plain fin-and-tube heat exchanger using genetic algorithm[END_REF] optimized four geometrical parameters in a plain fin-and-a tube heat exchanger using genetic algorithm. The objectives functions were the maximum total heat transfer rate and the minimum total pressure drop. The optimized heat exchanger has increased the heat transfer rate by 2.1-9.2% and the total pressure drop was reduced by 4.4 to 8 % for a Reynolds number ranging between 1200 and 14000.

Shape optimization

Bau [START_REF] Bau | Optimization of conduits' shape in micro heat exchangers[END_REF] optimized the shape of the cross section of the conduits in a micro heat exchanger by minimizing the maximum surface temperature of the conduit. He also found that the objective function could be further reduced by varying the width of the conduit in function of the axial coordinate. Hilbert et al. [START_REF] Hilbert | Multi-objective shape optimization of a heat exchanger using parallel genetic algorithms[END_REF] optimized the shape of blades in a tube bank heat exchanger, using a genetic algorithm suitable for multi-objective optimization. The two objectives functions are the temperature difference and the pressure difference at the optimization domain boundaries. The design variables are four geometrical parameters that define the shape of the blades. The optimization domain consists of four blades, whose position along the domain is fixed and all having the same shape. Gambit 2.1 was used for geometry and mesh generation, and Fluent 6.1 for solving the physical equations. The set of optimal solution formed the optimal Pareto front, which shows the conflict between heat transfer enhancement and pressure drop, as seen in Lee et al. [START_REF] Lee | Optimal shape and arrangement of staggered pins in the channel of a plate heat exchanger[END_REF] optimized the shape of pins and their arrangement in the channel in a plate heat exchanger. The design variables consisted of three parameters that define the shape of the pin and one parameter that defines the spacing between the pins. The optimization is based on the minimization of correlations for Nusselt number and friction factor, coupled in a single objective function. They used augmented Lagrange multiplier method to minimize the objective function. The optimum design variables defined the new shape and distribution of pins, which resulted in a 227.9 % enhancement in heat transfer and a 32.9% reduction in pressure drop, with respect to the baseline non-optimized design. Dlugosz [START_REF] Długosz | Parallel evolutionary algorithms in shape optimization of heat exchangers under thermomechanical loading[END_REF] used a sequential genetic algorithm for the optimization of a heat exchanger under thermomechanical load. He obtained three different shapes according to three optimization criteria: minimization of the maximal value of the temperature, minimization of the total volume of the heat exchanger and minimization of the maximal value of the equivalent stress. Figure 1.7: Pareto optimal front in shape optimization [START_REF] Hilbert | Multi-objective shape optimization of a heat exchanger using parallel genetic algorithms[END_REF] Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization [START_REF] Hilbert | Multi-objective shape optimization of a heat exchanger using parallel genetic algorithms[END_REF] Kwasi et al. [START_REF] Foli | Optimization of micro heat exchanger: CFD, analytical approach and multi-objective evolutionary algorithms[END_REF] used evolutionary algorithms to find the optimal shape of a separator between hot and cold side in a micro heat exchanger. The optimization was based on a multiobjective optimization technique, for simultaneously maximizing heat transfer and minimizing pressure drop. The design parameters that define the shape of the separator are represented by two "Non-Uniformal Rational B-splines, which consist of a number of points that define the polygon control shape". Optimization results were a set of optimal points that form the Pareto curve that shows the conflict between the objectives.

CFD-based optimization

Bougerard et al. [START_REF] Daniel | Optimisation de la géométrie de canaux d'échangeurs thermiques par simulation numérique[END_REF] used numerical simulation in size optimization of the canals of a heat exchanger. The study is performed under a fixed pressure difference equal to 40 Pa, for four types of canals: circular, square, isosceles right-angled triangle and equilateral triangle. The problem aims to find the optimal width of the heat exchanger canal using Gambit as the mesh generation software, Fluent as the CFD software and an optimization software, i-SIGHT-FD, in which multiple optimization algorithms are implemented. The optimization procedure consists of computing heat flux using CFD for a given value of the design parameter, then the optimization program generates a new value of the design parameter based on the value of the heat flux computed. Size optimization using CFD was also considered by Kwasi et al. [START_REF] Foli | Optimization of micro heat exchanger: CFD, analytical approach and multi-objective evolutionary algorithms[END_REF] in micro heat exchanger. The objective was to find the optimal aspect ratio, defined as the height of the channel to its width, for a constant and variable volume of microchannels.

Topology optimization

To the author knowledge, topology optimization was never applied in optimization of fluidto-fluid heat exchangers, beside a single research project [START_REF]Topology Optimization of Heat Exchangers[END_REF] where the method aims only to optimize the interfaces inside the optimization domain. On the other hand, topology optimization was widely applied in optimization of solid to fluid heat exchangers, as it will be seen later in this chapter. Next paragraph introduces in details the field of topology optimization, and the different methods used to solve this type of problem in order to select the best suitable strategy in applying topology optimization technique to the design of fluid to fluid heat exchangers.

Topology optimization 1.3.1. Introduction

Topology optimization was originally developed for the optimization of mechanical structures problems. The objective was to find the shape that uses minimum material while maintaining the mechanical stresses lower than an acceptable level (Figure 1.9). Topology optimization was defined by Bendsoe and Sigmund [START_REF] Bendsoe | Topology optimization: theory, methods and applications[END_REF] as shape optimization of continuum structures, which it should defines in every point in design space if there is a material in that point or not, the topology of the structure being not fixed a priori. Thus, from an initial field that is empty, full or in any intermediate state, the used control parameters allow the creations of holes and agglomerates of material without limitation, in order to find the best possible topology. The material distribution problem was the first application for topology optimization method in 1988 by Bendsøe and Kikuchi [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF]. Recently topology optimization concept was applied to a wide range of physical disciplines like fluids, heat transfer, acoustics, electromagnetic, optics. However the implementation of topology optimization methods is quite complex as it requires an algorithm capable of efficiently allocating and reallocating the material inside a predefined domain.

In topology optimization the domain space is discretized into finite meshes, also called cells, (cf. Figure 1.10) where each cell holds a design variable. The values of all design variables in all the cells define the shape of the entire structure. Thus the optimization problem is to find the optimal values of all the design variables, which minimize a certain objective function. Mainly there are two big families in topology optimization techniques: discrete approaches and continuous approaches, which stand respectively for methods in which the local parameter controlling the material in each cell has discrete values, or continuous bounded values. Those methods are gradient based optimization problems, thus the big number of design variables and the necessity to compute the derivative of the objective function with respect to the variables requires advanced mathematical techniques and software tools to ensure the convergence of the variables vector toward the optimal solution. Topology optimization has become a well developed area of research with many techniques developed to deal with numerical instabilities problems frequently encountered in topology optimization such checkerboards, meshdependencies and local minima. A represents an application of topology optimization on pure heat conduction problems. The final solution shows the optimal distribution of limited amount of high conductive material (in black) to evacuate the heat generated in low conductive material (in white) through a small patch (heat sink) located in the middle of the left boundary of the domain [START_REF] Marck | Topology optimization using the SIMP method for multiobjective conductive problems[END_REF]. Recently topology optimization was applied to fluid problems where the objective was to find the optimum layout of fluid channels to have the minimum pressure loss as shown in Figure 1.11.B [START_REF] Marck | Topology Optimization of Heat and Mass Transfer Problems: Laminar Flow[END_REF]. The white color represents the fluid that enters the domain from the entire left side and leaves it from a portion of right side, and the black color represents the solid. The shape of fluid channel is found by minimizing the total power dissipated by friction losses due to fluid flow through the domain. In [START_REF] Juan | Multi-objective optimization of a plain fin-and-tube heat exchanger using genetic algorithm[END_REF] a coupled thermal-fluid flow problem was solved using topology optimization as shown in Figure 1.11.C. It should be noted that in this type of problems the enhancement in heat transfer could be achieved only at the expense of pressure drop degradation and vice versa. Those examples make of topology optimization a promising methodology to solve optimization problems involving heat and mass transfer. A more detailed literature review on topology optimization of heat exchangers will be presented in the next paragraph.

Application in heat and mass transfer problems

Topology optimization methods

Problem formulation

Topology optimization problem is defined as follows: the domain is discretized into finite elements meshes; generally the same meshes used to solve the differential equations of the physical problem, in which a design parameter is stored. Those parameters are the design variables of the topology optimization problem. The design variables can take only discrete values; 0 or 1. In mechanical structure problem when the value of the design variable in a cell is 1, it means that the cell is made of material, while the 0 value means the cell is void. Similarly in fluid problems, 1 means the cell is solid and 0 the cell is fluid. The discrete nature of the problem makes it difficult to solve, therefore the design parameter is made continuous between 0 and 1 where the intermediate values of the design parameter represents a porous media or intermediate state.

The topology optimization methods aim to force the design parameter to take discrete values progressively, hence eliminate the grey regions and leads to a black and white domain where the design variables in each cell is either 0 or 1 (0-1 solution). However there also exist methods that can solve discrete combinatorial problems and are known as discrete approaches. We will consider a topology optimization problem where we're searching for the optimal distribution of two phases A and B. In discrete approaches the phase inside the cell is changed in a single step between A and B (hard-kill). In continuous approaches, a small quantity of phase A is replaced by phase B or vice versa in each step, until we reach at the end of the optimization process a cell completely made of A or B (soft-kill).

The general design problem is now defined as follows: find the material distribution; hence the optimal design variables, that minimizes an objective function , subject to a porosity constraint , and possibly other constraints . The optimization problem can be written in mathematical form as:

      0 0 1 min u( ), subject to z( ), 0, 1,...., 1 G ( ) 0 G u( ), 0, 1,...., 0 1, 1,...., l
N i i j i F G l K N jM iN                   (1.1)
where satisfies linear or non-linear state equations, represents the physical equations of the problem, the design variables, number of meshes and the maximum allowed porosity of material.

Topology optimization methods lie on three main parts: the direct solver of the physical problem (finite element, finite volume, …), the sensitivity analysis method (discrete adjoint, continuous adjoint, …) to compute the total derivative of the objective and constraint functions with respect to design variables, and a numerical optimizer to update the values of design variables at each iteration. Next we will present the most common methods present in the literature to solve problem (1.1).

Density method

The density method was initially developed by Bendsøe in 1989 [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF] in mechanical structure problem to find the optimum distribution of two materials, where one of them is void. The method uses continuous design variables taking values between 0 and 1. Those continuous variables are interpreted as material densities as they corresponds to void when having 0 value, to pure material when having 1 value and a composite material made of void and initial material when having intermediate values. The key part of this method is to introduce an interpolation function that computes various physical quantities in each cell, e.g. material stiffness in mechanical structure problems and thermal conductivity in heat conduction problems etc., as function of the continuous design variables:

  ( ) ( ) 1 ( ) i i A i B gg          (1.2)
Where is the density at node i, a physical quantity, and the physical quantities of the two phases at equal 0 and 1, and the interpolation function. The main challenge in density methods is the introduction of an interpolation function able to steers the solution to discrete 0-1 values and omits intermediate values, called a black-and-white design, with a real physical representations of composite materials corresponding to intermediate densities, known as grey material. A popular interpolation scheme to satisfy the above conditions is the penalized proportional "fictitious material" formula, known as the Solid Isotropic Material with Penalization (SIMP) defined as follows:

()

p ii g   (1.3)
where is the penalization parameter. Literature review shows that the optimal number that ensures good convergence to almost 0-1 solution is . However it is also preferred to begins the optimization with a low number of and increases it gradually. Bendsøe and Sigmund [START_REF] Bendsøe | Material interpolation schemes in topology optimization[END_REF] provided a physical justification of SIMP and showed that for physical realizability of intermediate densities elements is ensured. It is important to mention that without the presence of a volume constraint that limits material quantity the penalization effect in SIMP interpolation formula doesn't work out.

Stolpe and Svanberg [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF] introduced another interpolation function, Rational Approximation of Material Properties (RAMP) which is quite similar to SIMP in the principle of penalization effect to converge toward a black-and-white design. The interpolation function of RAMP is defined as follows [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF]:

  () 11 i i i g p      (1.4)
RAMP function was introduced after Stolpe and Svanberg demonstrated that SIMP interpolation scheme is non concave and the trajectory of global optimal solution in SIMP may be discontinuous [START_REF] Stolpe | On the trajectories of penalization methods for topology optimization[END_REF]. Pedersen [START_REF] Pedersen | Maximization of eigenvalues using topology optimization[END_REF] analyzed convergence problems associated with SIMP function at low density values, he suggested the use of density values higher than zero ( ) and low penalization parameters for densities lower than 0.1.

It should be noted that there is no general rule for the choice of penalization parameter that ensures a good convergence to 0 -1. Several other interpolation schemes have been developed for the same purpose of providing continuous interpolation with penalization effects (for example SINH method developed by Bruns [START_REF] Bruns | A reevaluation of the SIMP method with filtering and an alternative formulation for solid-void topology optimization[END_REF] using hyperbolic sinus function) and they are all stated as SIMP method.

Another approach used in density method is to add a concave function that serves as a penalty function that suppresses intermediate values in order to have a black-and-white-design.

The penalty function is added to the objective function or as an explicit constraint to the optimization problem [START_REF] Borrvall | Topology optimization using regularized intermediate density control[END_REF]. However this technique is not used in literature as much as the SIMP and RAMP interpolation functions.

Level set method

Level set methods are computational techniques introduced in 1988 by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] for tracking moving interfaces. The main idea of level set methods is to introduce a time and space dependent function that defines the interface between the two materials present in the topology optimization problem Г (e.g. fluid and solid, solid and void etc.) by the zero level contour of the called level set function. Having two material A and B in the optimization problem, the level set function is then defined as follows:

( , ) 0 if ( , ) 0 if ( , ) 0 if A B x t x x t x x t x                    (1.5)
where Ω A and Ω B represent respectively the domains of material A and material B. The evolution of level set function during the time defines the changes in the interface between the two materials. The level set function is updated from the solution of the Hamilton-Jacobi equation:

u. 0 t        (1.6)
where is the desired normal velocity on the boundary in which the zero level set propagates.

The main idea of level set as an optimization technique is to iteratively update the velocity , i.e. the direction, of the boundary propagation in order to decrease the objective function of the problem. Hence the geometry is defined through the optimization process by finding the optimal definition of material A -material B interface (e.g. solid -void interface in mechanical structure problems).

Level-set method is coupled in some applications with topological derivatives method to create new interfaces. Topological derivatives was first introduced by Eschenauer et al. [START_REF] Eschenauer | Bubble method for topology and shape optimization of structures[END_REF] for mechanical structure topology optimization and was recently used in fluid flow problems. The approach which is also known as 'bubble method' consists of finding the placement of insertion of infinitesimal hole in the design domain (small hole of phase A inside phase B sub-domain or vice versa). Coupling level set method with topological derivatives in fluid mechanics problem will be seen later in literature review.

Evolutionary approaches

Evolutionary approaches [START_REF] Huang | Evolutionary Topology Optimization of Continuum Structures: Methods and Applications Evolutionary Topology Optimization of Continuum Structures: Methods and Applications[END_REF] are considered among the most important branches of topology optimization using discrete design variables. Evolutionary Structural Optimization (ESO) was introduced by Xie and Steven [START_REF] Xie | A simple evolutionary procedure for structural optimization[END_REF] for optimization of mechanical structures. This method is based on the simple concept of removing inefficient material gradually from a structure until the constraint determining the volume of material in the design domain is no longer satisfied. Yang et al. developed the Bi-directional Evolutionary Structural Optimization (BESO), an improved version of ESO, which allows material to be removed and added simultaneously [START_REF] Young | 3D and multiple load case bidirectional evolutionary structural optimization (BESO)[END_REF]. The material removal and addition is based on the value of sensitivity number, where material is removed from elements with the lowest stresses (for mechanical structure problems for which the method was initially developed), and added into elements where the stresses are high. The sensitivity numbers in solid elements are estimated by the approximate variation of the objective function due to the removal of elements whereas in void elements sensitivity numbers are set to zero.

The same approach was used in topology optimization of conductive heat transfer problems, allowing adding or removing elements in function of their effective contribution to the heat transfer [START_REF] Li | Evolutionary topology optimization for temperature reduction of heat conducting fields[END_REF], [START_REF] Li | Shape and topology design for heat conduction by Evolutionary Structural Optimization[END_REF]. In comparison with SIMP method, BESO could be categorized as a discrete version of SIMP scheme. Furthermore Huang and Xie introduced a version of BESO [START_REF] Huang | Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials[END_REF] called soft-kill BESO in which they use a power law interpolation parameterization (SIMP) for the computation of gradients to update the design variables. In soft-kill BESO the design variables still have discrete values only but the gradients are derived from continuous variables in function of penalization parameter . Hence, the original BESO method, called hard-kill BESO in which the gradient have discrete values similar to design variables, is a special case of the softkill BESO method where the penalty exponent tends toward infinity.

Many works in topology optimization of mechanical structure problems used the evolutionary approaches with many techniques developed to deal with the difficulties encountered due to the discrete aspect of the problem. In pure heat conduction, evolutionary approaches were also applicable but the results showed that the method leads to local optimum. Therefore evolutionary approaches have not gained interest in topology optimization of fluid mechanics.

Topology optimization in heat and mass transfer problems, case of two fluids

As already mentioned, the aim of this thesis is to apply topology optimization technique on design and optimization of fluid to fluid heat exchangers. This requires three main steps:

1. Formulation of a suitable topology optimization method able to consider the distribution of three phases in the domain, two fluids and one solid. 2. Separation of fluids sub-domains to allow each fluid to connect its predefined inlet sections to outlet sections while avoiding any fluid mixture inside the domain. Next a literature review on each step of the thesis objective will be presented, in order to select the best strategy regarding every step. It should be noted that to the author knowledge, no work have been done in phase separation, hence literature review is limited on the first and third steps. The application of topology optimization was first performed by Borrvall and Petersson [66] on Stoke flows. Stokes flow is a type of fluid flow where the advective inertial forces are insignificant compared to viscous forces, thus having low Reynolds number (Re << 1). Borrvall and Petersson provided mathematical proofs of existence of an optimal solution for topology optimization in fluid using density approach and an appropriate penalized interpolation function. The optimization problem consists of minimizing the dissipated power in a fluid domain; the total volume of fluid should not exceed a maximum value considered as a constraint. Figure 1.12.A shows the optimal solution of a pipe bend which corresponds to a straight line having the minimum distance between inlet and outlet [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF]. Figure 1.12.B and Figure 1.12.C show that when the domain is square the optimal shape of pipes is two straight parallel lines, whereas when the domain is rectangular with length 1.5 times bigger than the width, the minimum pressure drops is reached by joining the two pipes in the middle [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF].

Literature review on topology optimization in mass transfer problems

Aage et al. [START_REF] Aage | Topology optimization of large scale stokes flow problems[END_REF] solved the same problem for large scale domains using parallel computations techniques. Gresborh-Hansen et al. [START_REF] Gersborg-Hansen | Topology Optimization of Mass Distribution Problems in Stokes Flow[END_REF] used outflow rate as optimization target.

Guest et al. [START_REF] Guest | Topology optimization of creeping fluid flows using a Darcy-Stokes finite element[END_REF] proposed to treat the solid phase in the optimization domain as a porous medium with flow governed by Darcy's law. Fluid phase and solid phase are then treated in a single equation created by combining Stokes and Darcy equations. The resulted Darcy-Stokes equation leads to an appropriate modeling for no-slip boundary condition along the fluid-solid interface inside the domain.

Wiker et al [START_REF] Wiker | Topology optimization of regions of Darcy and Stokes flow[END_REF] also treated Darcy-Stokes topology optimization problem. Their work included solution to an area to point drainage problem, where the goal was to transport all fluid out from the domain through a single part of the boundary with minimal power consumption possible. They also studied the impact of many geometrical and mathematical parameters on the final solution. Contrary to [START_REF] Guest | Topology optimization of creeping fluid flows using a Darcy-Stokes finite element[END_REF], their results showed that regularization is needed in topology optimization of fluid problems to avoid numerical problems.

Evgrafov et al. [START_REF] Evgrafov | The Limits of Porous Materials in the Topology Optimization of Stokes Flows[END_REF] stated that the problem of Darcy-Stokes flow in topology optimization admit solutions even if the limiting zero and infinite permeabilities of Darcy's law are included in the design domain.

Topology optimization of fluids in Stokes flow in Borrvall and Petersson work was extended by Gersborg-Hansen et al. [START_REF] Gersborg-Hansen | Topology optimization of channel flow problems[END_REF] to Navier-Stokes flows by describing a topology optimization method for steady, incompressible laminar viscous flows at low to moderate Reynolds number with inclusion of inertia effects which made the flow problem nonlinear. The use of high-level programming frameworks for topology optimization of steady-state Navier-Stokes flow problems were considered in [START_REF] Olesen | A high-level programming-language implementation of topology optimization applied to steady-state Navier-Stokes flow[END_REF].

The articles on topology optimization for Navier-Stokes flow stated so far refers only to steady state problems. However unsteady flow could also be implemented in topology optimization methods where the feasibility of the problem is well demonstrated in [START_REF] Deng | Topology optimization of unsteady incompressible Navier-Stokes flows[END_REF]. The results showed that the final optimum design of unsteady Navier-Stokes flow is influenced by the dynamic effect and Reynolds number. Kreissl et al. [START_REF] Kreissl | Topology optimization for unsteady flow[END_REF] studied the feasibility of density approach for optimizing the topology of unsteady flows. Results showed that structures optimized for unsteady flow differ increasingly from corresponding steady-state design as the problem becomes more unsteady.

Deng et al. [START_REF] Deng | Topology optimization of steady and unsteady incompressible Navier-Stokes flows driven by body forces[END_REF] added to the unsteady problem flows driven by body forces. The physical body forces in the equations were penalized using power law approach, and results showed that for both steady and unsteady flows, optimal topology designs depend on the type of the body forces.

Level-set

Zhou and Li [START_REF] Zhou | A variational level set method for the topology optimization of steadystate Navier-Stokes flow[END_REF] performed 2D and 3D numerical experiments in topological design using level set method, and proved a relatively good agreement of their results with those obtained by density methods. Figure 1.13 and Figure 1.14 show respectively the topology optimization results using level set method of bend pipe and a diffuser [START_REF] Duan | Shape-topology optimization for Navier-Stokes problem using variational level set method[END_REF], which was previously studied using density method [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF], [START_REF] Gersborg-Hansen | Topology optimization of channel flow problems[END_REF]. However, even if for some specific examples level set method showed agreement with density method, the performance of this technique is still significantly limited by its incapacity of creating new holes in the design domain. Hence the topology optimization procedure is strongly dependent on the initial guess. However this issue was first solved in topology optimization of mechanical structure problems by incorporating topological derivative method into level set method to reduce its dependency on the initial guess [START_REF] Yu | Topology optimization with Level Set Method Incorporating Topological Derivative[END_REF], then it was applied to fluid problems by various authors [START_REF] Deng | Topology optimization of steady Navier-Stokes flow with body force[END_REF], [START_REF] Challis | Level set topology optimization of fluids in Stokes flow[END_REF].Topological derivative indicates then the location where new holes should be nucleated using two different strategies: even including topological sensitivity information in the evolution equation of level set method, or using this information at discrete places in the optimization algorithm [START_REF] Challis | Level set topology optimization of fluids in Stokes flow[END_REF].

Literature review on topology optimization in heat and mass transfer problems

Topology optimization of convective heat transfer was first solved by Bruns in 2007 [START_REF] Bruns | Topology optimization of convection-dominated, steady-state heat transfer problems[END_REF]. The author introduced a method for solving conductive and convective non linear steady state heat transfer problems by adding a convection term to heat transfer equation but without evaluating the fluid motion.

Matsumori & al. [START_REF] Matsumori | Topology optimization for fluidthermal interaction problems under constant input power[END_REF] used the density approach to solve fluid thermal interaction topology optimization problems. Instead of taking two objective functions for power dissipation and heat transfer in a multi-objective optimization approach, they proposed a method to keep the input power constant within the optimization process and solved the optimization problem by maximizing only heat transfer objective function. The optimization problem was formulated for two types of heat conditions: temperature dependent heat sources and temperature independent heat sources. Okkels et al. [START_REF] Okkels | Applications of topology optimization in the design of micro-and nanofluidic systems[END_REF] optimized micro-cooling system using a high level programming language for the implementation of non linear topology optimization using Matlab and the commercial finite element software Femlab.

Dede [START_REF] Dede | Multiphysics topology optimization of heat transfer and fluid flow systems[END_REF] solved heat and mass transfer problem in topology optimization problem by coupling. the mean temperature and power dissipation simultaneously using a weighted sum approach for multi-objective optimization. The method was applied to the design of hierarchical microchannel system [START_REF] Dede | Experimental Investigation of the Thermal Performance of a Manifold Hierarchical Microchannel Cold Plate[END_REF] and to a submerged jet impingements to remove heat from an electronic package that generate heat to ensure that the maximum device temperature does not exceed an allowable limit [START_REF] Dede | Multiphysics optimization, synthesis, and application of jet impingement target surfaces[END_REF].

In [START_REF] Dede | Optimization and Design of a Multipass Branching Microchannel Heat Sink for Electronics Cooling[END_REF] Dede used two different interpolation schemes, one for the calculation of thermal conductivity and another one for the calculation of the inverse permeability using the same design variable. He applied the same procedure as in [START_REF] Dede | Multiphysics topology optimization of heat transfer and fluid flow systems[END_REF] for optimizing microchannels for the cooling high heat generation electronics devices. Yoon [START_REF] Yoon | Topological design of heat dissipating structure with forced convective heat transfer[END_REF] studied the minimization of thermal compliance in a heat dissipating structure under constant mass flow using SIMP formula for the interpolation of material properties that varies between solid and fluid phases. The method was then extended and applied to the design of an electro-thermal-compliant actuator device [START_REF] Yoon | Topological layout design of electro-fluid-thermal-compliant actuator[END_REF].

Koga et al. [START_REF] Koga | Development of heat sink device by using topology optimization[END_REF] developed a heat sink device for small scale applications considering low velocities and low Reynolds number neglecting inertial forces, thus the study is limited to Stokes flow. The results that combined multi-physics objective function involving the pressure drop and heat transfer performance showed that when the weight of heat transfer criterion is increased in the multi-objective function, small auxiliary channels appear inside the domain. These ramifications increase the heat distribution which leads to a more efficient heat transfer over the domain but increases the fluid flow pressure drops. Yaji et al. [START_REF] Yaji | A topology optimization method for a coupled thermal-fluid problem using level set boundary expressions[END_REF] developed a method for solving the coupled thermal fluid topology optimization problem using level set method in contradiction with the methods stated above in which they used material interpolation schemes.

As mentioned previously, topology optimization are in generally gradient-based optimization methods, because of the high number of design variables This is seen from literature review, in which all works presented used gradient-based numerical optimization algorithms, especially the method of moving asymptotes. However, Yoshimura [START_REF] Yoshimura | Topology optimization of fluid problems using genetic algorithm assisted by the Kriging model[END_REF] used genetic algorithm instead of gradient based optimization algorithm to solve topology optimization problems in fluid flow using level set approach. To reduce the computational cost caused by the evaluation of the objective function, he assisted the genetic algorithm by Kriging surrogate model to estimate the objective function values and their uncertainties. However, due to the low number of design variables, the results were less accurate than those obtained in conventional methods that use gradient information in the optimization process. 1.1 provides a summary on applications of topology optimization on heat and mass transfer in mono-fluid domain.

Article

Conclusion

After being initially developed for mechanical structure optimization problems, topology optimization approaches are now widely used to various flow optimization problems. Optimization for flow problems was initially limited to low Reynolds numbers and steady state (Stokes flow) without taking into consideration inertia effects. Then various authors have extended the optimization procedure to cover wider range of Reynolds number, inertia effects (Darcy-Stokes and Navier-Stokes flows), non-uniform body forces and unsteady state flows. Density approach, level set, topological derivatives and evolutionary approaches especially BESO method are the most used approaches. The main difference between these different approaches is the design variables update schemes. Evolutionary approaches use discrete design variable whereas the rest of methods uses bounded continuous design variables. It's important to point that filtering and smoothing techniques are applied in the majority of topology optimization approaches to deal with numerical problems and instabilities such as convergence toward local minimum, mesh dependencies solutions and checkerboard problems.

Despite gaining attention in mechanical structure and pure diffusion heat transfer problems, evolutionary techniques have not been considered in fluid flow problems according to literature review. Level set are found to be appealing for flow problems due to its results in 2D and 3D numerical simulations for various type of flows. However, level set methods can only evolve from existing boundaries and is not able to generate new holes, which means unable to generate new channels in fluid topology optimization. This is considered as a conceptual drawback for the method, especially if it will be used for topology optimization of heat exchangers. Nucleation of new holes in level set method was possible by combining level set method with topological sensitivity information method. This combined method was applied and tested by various authors. Results show that the combined method still strongly dependent on the starting guess [START_REF] Sigmund | Topology optimization approaches[END_REF].

The density method is completely independent from the initial guess and the generation of channels and complex structures is only dependent from the objective and constraint functions. Moreover Table 1.1 shows that density method was applied on majority of problems that coupled heat and mass transfer. Despite requiring high computation time, density method appeared to be the most suitable method to extend the application of topology optimization in fluid mechanics to bi-fluid domain, which to the author knowledge was not considered before

Outline of the research

As concluded in last paragraph, the density approach is the most suitable method for topological optimization of bi-fluid heat exchangers and it will be considered in this work. The general algorithm of density method we applied to heat and mass transfer problems is composed of three main steps:

 The CFD solver using a finite volume method.

 The sensitivity analysis based on discrete adjoint method.  The optimizer that computes the new values of the design variables on the basis of objective function, the constraints functions and their derivatives. The method of moving asymptotes was chosen.

The remaining of the document is divided as follows:

 The detailed algorithm of optimization method presented above, and detailed development of each part of the method are presented in chapter 2. Two different formulations will be compared, one using a single design variable in each design cell and the second one uses two design variables in each design cell which double the number of variables of the problem.  Fluid separation will be considered in chapter 3. Three different methods are formulated and implemented in the optimization algorithm.  In chapter 4 maximization of heat transfer between the two separated fluids will be considered.  Finally the conclusion and perspectives in chapter 5.

Chapitre 2: Optimisation topologique bi-fluide 1. Introduction

Ce chapitre introduit une méthode d'optimisation topologique basée sur l'approche de la densité appliquée aux problèmes de mécanique des fluides pour l'écoulement laminaire bifluidique. L'algorithme général de la méthode sera présenté avec des explications détaillées sur la mise en oeuvre de chaque étape de la méthode. Deux formules d'interpolation reliant les variables de conception et les variables physiques des deux fluides et du solide sont générées, testées et comparées afin de sélectionner la plus appropriée pour la résolution du problème bifluide. Les fonctions de contraintes qui limitent le volume maximum de fluide seront également décrites.

Formulation du problème mathématique

Le problème d'optimisation est maintenant défini comme suit: trouver les valeurs optimales des variables de conception qui définissent la distribution optimale des trois phases, fluide solide 1 et fluide 2, afin de minimiser la dissipation de puissance du fluide dans le domaine, sans dépasser les porosités maximales autorisées de chaque fluide. La fonction objectif dépend explicitement des champs de vitesses et de pression. Le calcul des paramètres physiques dans les mailles de conception utilisant la fonction d'interpolation pénalisée, fait que la fonction objectif dépend implicitement de la variable de conception η. Le problème d'optimisation est résumé dans (2.8).

Modélisation de l'écoulement et de transfert thermique

Les équations physiques qui représentent l'écoulement du fluide et le transfert thermique par conduction et transport, sont l'équation de continuité (2.4), de Navier-Stokes (2.5) et l'équation d'énergie (2.7).

Pour modéliser les mailles solides et fluides en utilisant la même équation de quantité de mouvement, Borrvall et Petterson [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] ont proposé d'ajouter à cette équation un coefficient inverse de perméabilité α similaire à la loi de Darcy dans un milieu poreux. Ainsi pour une maille fluide α sera égal à 0, et une pour une maille solide α sera égal à une valeur suffisamment grande pour le solide, donc à la limite, la vitesse dans les régions solides sera nulle. Par conséquent, lorsque la variable de conception prend une valeur intermédiaire, la perméabilité inverse prendra une valeur entre zéro et sa limite maximale, de sorte que la cellule de conception correspondante sera un milieu poreux dans lequel la vitesse du flux est ralentie en fonction de la quantité de solide dans la cellule poreuse. Le terme de perméabilité mathématiquement inverse agit comme un terme d'absorption de vitesse. Physiquement, cela pourrait être vu comme une force interne du corps, ou une force de friction proportionnelle à la vitesse du fluide entre le fluide et un petit obstacle solide à l'intérieur du flux Il convient de noter qu'il n'est pas nécessaire d'insérer le terme de diffusion pour la modélisation de l'écoulement de plusieurs fluides à travers un milieu poreux. En fait, dans la méthode d'approche de densité appliquée aux problèmes d'écoulement de fluide, dans un état intermédiaire, le fluide 1 et le fluide 2 formeront un troisième fluide équivalent, qui forme à son tour avec le solide un matériau poreux unique. Les propriétés physiques sont ainsi interpolées entre le solide, le fluide 1 et le fluide 2. Ainsi, la cellule entière a un seul vecteur de vitesse et il n'y a pas de vitesse de circulation pour chaque fluide.

Algorithme d'optimisation

L'approche de densité consiste à faire varier de façon continue le paramètre de pénalisation p de la fonction d'interpolation, jusqu'à ce qu'il atteigne sa valeur la plus appropriée pour assurer une convergence en douceur vers la solution en noir et blanc. L'algorithme d'optimisation est donc divisé en deux boucles:  Une boucle principale (Figure 2 En optimisation topologique, plusieurs problèmes numériques peuvent être rencontrés tels l'apparition d'une succession de régions solides et fluides dans des mailles adjacentes. Ce problème est résolu en utilisant les techniques de filtrage numériques. Dans cette thèse, des filtres de densité et de sensibilité sont utilisés. Le filtre de densité (équation (2.53)) consiste à remplacer la variable de conception dans chaque maille du domaine par une moyenne pondérée en fonction des variables de conception dans les mailles voisines. Ce filtre est appliqué entre la 5 ème étape et la 1 ère étape de boucle interne. D'autre part, le filtre de sensibilité (équation (2.55)) est la redéfinition de la dérivée par rapport à la variable de design dans une maille en fonction des valeurs de la dérivée dans les mailles voisines. Ce filtre est appliqué sur le champ de dérivé entre la 4 ème et 5 ème étape de la boucle interne. Le voisinage d'une maille est défini en général comme les éléments dont les centres sont à une distance prédéfinie du centre de la maille concernée (équation (2.54)).

Fonction d'interpolation pénalisée a) Fonction d'interpolation mono-eta

Le premier schéma d'interpolation utilise une seule variable de conception dans chaque maille, de même que dans les problèmes mono-fluide. Ce schéma est basé sur la fonction de distribution normale (2.9). Dans cette fonction est la variable de design, est l'abscisse de la pique de la courbe et le paramètre qui détermine la forme de la courbe (Figure 2.4). Yin et al [START_REF] Yin | Topology optimization of compliant mechanisms with multiple materials using a peak function material interpolation scheme[END_REF] ont utilisé la même fonction de distribution dans l'optimisation topologique des problèmes de structure multimatériaux. Ainsi la fonction d'interpolation sera la superposition de trois fonctions normales, chacune correspondante à l'une des trois phases (fluide 1, fluide 2 ou le solide), et chacune multipliée par la propriété physique de la phase qui lui correspond. Il reste à déterminer la distribution des valeurs discrètes de correspondant aux phases fluide et solide, c.à.d. déterminer l'abscisse de la pique de chaque phase. La distribution qui assure la représentation physique la plus pratique des valeurs intermédiaires de

Résultats a) Formulation mono-eta

La configuration présentée dans la Figure 2.12 est un domaine carré de 36 mm de longueur ayant une entrée et une sortie pour chaque fluide. La porosité maximale de chaque fluide est égale , ce qui est équivalent au volume de fluide requis par un canal droit reliant l'entrée et la sortie de chaque fluide. La masse volumique, la viscosité dynamique et le coefficient de perméabilité sont égaux pour les deux fluides. Afin d'analyser le comportement de la formulation mono eta en ce qui concerne le point stationnaire, quatre études ont été réalisées, chacune correspondante à une initialisation différente de la variable de conception. La Figure 2.13 montre les solutions finales de chaque cas d'étude. Dans le cas où , aucun changement n'a eu lieu dans les valeurs des paramètres de conception en raison de la dérivée nulle pour , ainsi la dérivée de la fonction objectif sera nulle. Pour , l'algorithme distribue toute la quantité permise du fluide 2, représentée par la couleur rouge, entre les sections d'entrées et sorties des deux fluides. Cependant, les deux tuyaux formés par le fluide 2 ne remplissent pas la section complète entre les sections d'entrée et de sortie de chaque tuyau, en raison de la contrainte de porosité du fluide 2 qui est égale seulement au volume d'un tuyau complet. La quantité de fluide 1 représentée par la couleur bleue, qui a été ajoutée aux étapes ultérieures du processus d'optimisation n'a pas pu être ajoutée à l'endroit optimal dans les régions d'écoulement afin de minimiser la perte de charge. Cela montre la difficulté d'ajouter du fluide 1 lorsque la valeur initiale de la variable de conception est supérieure à 0,5, à cause de la forme de la courbe de la fonction d'interpolation. Pour , le même comportement est observé mais cette fois pour une convergence optimale vers le fluide 1. Enfin pour le cas où dans les régions au voisinage de l'entrée et la sortie du fluide 2 et dans les régions au voisinage de l'entrée et la sortie du fluide 1, deux tuyaux sont formés, chacun comprenant la quantité maximale totale permise de l'un des deux fluides. Cela montre que la formulation mono-eta est fortement dépendante sur l'estimation initiale, comme déjà discuté ci-dessus.

b) Formulation bi-eta

La même configuration que le cas précédent est reconsidérée, cette fois le problème est résolu par la formulation bi-eta. La solution finale dans la This chapter introduces a topology optimization method based on density approach applied on mechanical fluid problems for bi-fluidic laminar flow. The general algorithm of the method will be presented with detailed explanation on the implementation of the direct solver based on finite volume method and sensitivity analysis using discrete adjoint approach. Two interpolation formulas relating the design variables and the physical variables of the two fluids and the solid are generated, tested and compared in order to select the most suitable one for the problem. Constraints functions that limit the maximum fluid volume will be also described.

Problem formulation

The domain is made of three subdomains: solid subdomain , first fluid subdomain representing the cold stream , and the second fluid subdomain representing the hot stream . The subdomains represented in Figure 2.1 verify the materials conservation equation:

. It should be noticed that in the final solution there is no intersection between the three subdomains. Boundary conditions could be: Fluid is assumed to be Newtonian and incompressible under steady state laminar regime. Furthermore the whole domain can be subject to a constant heat generation rate.


   0 T n    2 f  1 f  Ω s Ω s Ω s

Fluid flow modeling

The fluid flow problem is solved using the Navier-Stokes equations, described by the continuity equation (2.1) and momentum equation (2.2) as follows [START_REF] White | Fluid Mechanics[END_REF]:

  . u 0 m m t        (2.1) 2 u u. u u m PY t              (2.2)
where is the density of fluid mixture supposed constant (in intermediate state both fluids can exist in the same cell), and are respectively velocity and pressure fields, the fluid dynamic viscosity and the fluid body forces.

For modeling the solid regions and fluid regions using the same momentum equation, Borrvall and Petersson [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] proposed to add to the momentum equation an inverse permeability coefficient similar to the Darcy's law within a porous media. Inverse permeability coefficient ranges from zero value for fluid to a sufficiently large value for solid, thus at the limit, velocity in solid regions goes to zero. Hence when the design variable takes intermediate value, the inverse permeability will take a value between zero and its maximum limit, so the corresponding design cell will be a porous media in which the velocity of the flow is slowed down in function of the quantity of solid in the porous cell. Mathematically inverse permeability term acts like a velocity absorption term. Physically this could be seen as internal body force, or a friction force proportional to fluid velocity between the fluid and a small solid obstacle inside the flow. in equation (2.2) is expressed as follows:

u Y   (2.3)
Finally considering a steady state flow the derivative of the velocity with respect to time in equation (2.2) will be equal to zero, and considering an incompressible Newtonian fluid the derivative of density with respect to time in equation (2.1) is also null. Combining equation (2.3) with equation (2.1), the fluid flow equations are expressed as follows:

.u 0

 (2.4)   2 ( ) u. u ( )u ( ) u P             (2.5)
It should be noted that fluid density and dynamic viscosity values in equations (2.4) and (2.5) are function of the properties of the first and second fluids (equivalent fluid is mixture of both fluids). Inverse permeability is function of equivalent fluid properties, solid material properties or intermediate porous state between fluid and solid. The variables , and are then function of the variable design according to density approach.

Theoretically in order to have zero velocity in solid regions, should be for solid and for fluid to model pure fluid flow without obstacles. In practice the two limits cannot be applied, so should be bounded so that the lower value is and upper value is , as suggested in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF]. The inverse permeability value to enforce zero velocity in solid regions depends on Reynolds number. Kreissl et al. [START_REF] Kreissl | Topology optimization for unsteady flow[END_REF] studied the variation of velocity with respect to for different Reynolds number and showed that when Reynolds number decreases larger value of are needed to force the velocity to go to zero in solid regions. However results showed that is sufficient to ensure zero velocity in a solid cell for all Reynolds number.

It should be noted that it isn't necessary to insert the diffusion term in left side of continuity equation (2.1) similarly to multi-component fluids flow through a porous medium. Actually in density approach method applied to fluid flow problems, in intermediate state both fluid 1 and fluid 2 will form a third equivalent fluid, which in his turn form with the solid a single porous material. The physical properties of the porous are interpolated between the solid, fluid 1 and fluid 2. Hence the entire cell has a single velocity vector and there isn't a flowing velocity for each fluid apart.

Heat transfer modeling

For an isotropic porous medium where the radiative effects, viscous dissipation, and the work done by pressure changes are negligible, assuming a thermal equilibrium between the solid and fluid phases, the energy equation is expressed as follows [START_REF] Nield | Convection in Porous Media[END_REF]:

      u. . P P m m m f T C C T k T q t          (2.6)
where is the density, is the specific heat, is the thermal conductivity, the heat production per unit volume, the temperature and the time. Subscript corresponds to a mixture of solid and two fluids. Subscript corresponds to the mixture of two fluids without solid. In real physical application and are weighted mean values of solid and fluid properties function of porosity.

The evaluation of the (equivalent) thermal conductivity of the porous medium is a more complicated task and depends in a complex fashion on the geometry of the medium [START_REF] Nield | Convection in Porous Media[END_REF]. For the simple case where the medium is homogeneous, the overall thermal conductivity is evaluated by the weighted average value of the solid and fluids conductivities. This is considered as the upper bound of the actual overall thermal conductivity according to [START_REF] Nield | Convection in Porous Media[END_REF].

In the present work, porous medium exists for intermediate design only and can then be considered as a homogeneous material. Therefore homogeneous porous material properties will be calculated using the penalized interpolation function instead of weighted average value function. All the physical quantities of equation (3.6) are evaluated in the same manner. It should finally be noted that in steady state condition the term will be equal to zero. Equation (3.6) becomes:

    ( ) ( ) u. . ( ) ( ) P C T k T q           (2.7)

Optimization problem

The topology optimization problem aims finding the optimum distribution of , and subdomains defined by a vector of local design parameter , in order to maximize or minimize an objective function The objective function depends explicitly on u, and respectively the flow velocity, pressure and temperature fields computed by equations (2.4), (2.5) and (2.7). depends implicitly on through the physical quantities computed by the interpolation function, which are present in Navier-Stokes and energy equations, and also could depend explicitly on the design parameter . Finally the problem is subjected to two porosity constraints for fluid 1 and fluid 2 respectively. The optimization problem (1.1) is now defined as follows:

  11 22 min z( ), with z( ) (u, , ) subject to Equations (2.4) (2.5) (2.7) ( ) ( ) ff ff F f P T G G          (2.8)
where and represent respectively the maximum allowed porosity of fluid 1 and fluid 2 in the design domain.

Algorithmic scheme

Topology optimization method consists of modifying the three sub-domains in order to minimize the objective function in problem statement (2.8). The density approach discussed in detail in paragraph 1.4.2 consists of changing continuously a local design parameter stored in each element over the discritized domain. The impact of each element subdomain's belonging on the objective function is provided by the objective function's gradient and serves as the convenient information to the numerical optimization method to converge design parameter vector towards its optimal solution. It should be noted that density approach method consists of varying continuously the penalization parameter of the interpolation function, until it reaches its most convenient value that ensures smooth convergence to black-and-white solution. The different steps required to apply density approach using the penalization effect is divided into two loops:

 Main loop (cf. Figure 2.
2) consists of repeating the inner loop by changing the penalization parameter progressively. The number of iterations of the main loop is evaluated by means of an integer ranging from 0 to a prefixed number allowing penalization parameter to reach successively its most convenient value. is the value of the objective function computed in each iteration of the inner loop. and are the residuals of the objective function in inner and outer loops respectively, to test the convergence of the problem. and are initialized with , which stands for a number larger than the objective function, in order to prevent the convergence from the first iteration. This loop is repeated until reaching the global convergence. 1. Computing porous media properties using the interpolation function and the design parameters from the previous iteration. 2. Solving the direct problem using finite volume method to find velocity, pressure and temperature fields. 3. Evaluation of the objective function. 4. Evaluations of the derivative of the objective function with respect to the design parameter using discrete adjoint approach. 5. Numerical optimization algorithm: Create the subproblem of the method of moving asymptotes, which consists of finding an approximate convex function for the objective and constraint functions, then solving the subproblem using the interior point method..

Several regularization functions (filters) are used at different level of the algorithm. These regularization functions allow better convergence of the optimization problem.

Interpolation functions with penalization

Topology optimization method applied to fluid flow problems have to decide whether to place solid or fluid in each element of the domain. Fluid and solid phases are directly defined using a discrete design variable stored in the design cell. The density approach with penalization aims to replace the discrete variables by continuous variables. Penalization function coupled with porosity constraints will then forces these variables to converge to desired discrete values.

The physical properties of each element of the porous media as well as its permeability are defined as a function of this continuous design variable using an interpolation function.

In mono fluid optimization problem the design variable varies between 0 and 1, where 0 corresponds to solid phase, 1 corresponds to fluid phase and intermediate values of design variable corresponds to a porous medium for the case in which fluid and solid exist simultaneously in the same design cell.

The actual problem needs an interpolation function able to assign physical quantity in each cell in function of the properties of the three phases : solid phase, fluid 1 phase and fluid 2 phase. In this work, two interpolation schemes were developed and tested: The peak function based on normal distribution function using a single design variable in each cell, and the multi-phases interpolation function that uses two designs variables in each cell.

CFD Solver

Mono-eta interpolation function

The first interpolation scheme uses a single design variable in each cell, similarly to two materials phases case (fluid and solid). The interpolation scheme is based on the normal distribution function given as follows:

  2 2 exp 2 N p        (2.9)
where is the design variable, the location parameter of the peak and is a scale parameter. and determine the shape of the normal distribution function, as shown in Figure 2.4. Yin et al [START_REF] Yin | Topology optimization of compliant mechanisms with multiple materials using a peak function material interpolation scheme[END_REF] used the same function for topology optimization of compliant mechanisms with multiple material. The constitutive linear material property tensor is the sum of the properties of materials from which the mechanism is formed multiplied by the normal distribution function, as follows:

  2 2 1 exp 2 n im im m m p           (2.10)
When the design variable in a particular cell approaches , which correspond to a phase , its normal distribution function approaches 1, which means the physical property in this cell will be equal the physical property of 's phase. In case of interpolation of 2 fluids and 1 solid, considering a design cell in the domain, : hence: and . This distribution of each phase's peak location along axis was chosen to have the most convenient physical representation of intermediate values of : for we will have a porous media solid -fluid 1 and for we will have porous media solid-fluid2. Finally the mono-fluid interpolation scheme is expressed as follows:

    22 2 12 2 2 2 12 0.5 1 exp exp exp 2 2 2 ii i i f s f f s f p p p                                  (2.12)
Equation (2.12) is drawn in

  2 1 1 2 1 exp 2 f f c Gd Np             (2.13)   2 2 2 2 1 1 exp 2 f f c Gd Np              (2.14) 0 0.2 0.4 0.6 0.8 1 Design variable η 1 f  2 f  s  0.05 p  0 0.2 0.4 0.6 0.8 1
Design variable η s  1 2 , f f   0.05 p 
where and are respectively the maximum allowable volume of fluid 1 and fluid 2 in the optimization domain. It should be noted that is the number of design variables , and being respectively the number of horizontal and vertical design cells for a 2D design domain.

A special attention should be given to the choice of the penalization parameters of the interpolation function and the constraint functions. Figure 2.6 shows the constraint functions for both fluid 1 and fluid 2 for Both functions goes approximately to 0 at which corresponds to the peak point of solid phase, thus this value of penalization parameter could be considered as the ideal value for a good representation of fluid 1 and fluid 2 volume percentage in a cell in function of design variable . Penalization parameters of interpolation function (2.12) and which correspond respectively to solid, fluid 1 and fluid 2 peaks could be taken equal to a single value . In order to ensure the convergence of design variables in all cells to 0, 0.5 or 1 and to omit all intermediate values of , the interpolation function of inverse permeability coefficient should have a continuous negative slope between and and between and , thus the minimum value of should strictly be at and . The slope of the interpolation function should also not decrease dramatically near 0.5, to prevent reaching the minimum of permeability at or , otherwise high level of intermediate material will remain in the final solution. This requires a large value of as shown in Figure 2.4.

Another important feature of the interpolation function that should be discussed is the stationary point at that have null derivative . The entire optimization algorithm that search for the optimal value of depends on the slope of the local point on the interpolation function. Thus the zero slope of represents a major difficulty for the algorithm to cross the stationary point and make transition from fluid 1 phase to fluid 2 and vice versa. Hence the algorithm will be strongly dependent on the initial guess, for example if the initial value of is taken somewhere between 0 and 0.49 the algorithm will have difficulty to add fluid 2 in optimal location, similarly if initial is taken between 0.51 and 1 the algorithm will easily add fluid 2 and will have difficulty to add fluid 1. This behavior will be clearly shown later in application part where weak and tough points of mono-eta interpolation function will be discussed.

Bi-eta interpolation function 2.4.2.1. Mono fluid interpolation function

In [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] Borrvall and Peterson introduced the first interpolation function to fluid flow problems, they also provided the mathematical proofs that function allow the convergence to an optimal solution. Many authors have latter adopted this interpolation scheme. According to this function, the physical properties in each element of the domain vary continuously as a function of a design variable ranging in interval [0,1]. The interpolation convex function is defined as follows:

  1 ( , )
with 0 1 and 0

s f s p p p p                 (2.15)
where is a scalar standing for physical properties (for example inverse permeability coefficient) with and the corresponding values for solid and fluid phases respectively. The parameter is a penalty parameter added to adjust the convexity of this function as presented in Figure 2.7, and plays a decisive role in the convergence of the entire method. When the penalization parameter tends to zero, the penalization function become more convex and the optimization algorithm allows higher number of gray cells. At high values of penalization parameter (

), the interpolation function will have pseudo linear shape and will reduce the number of gray cells. The level of grey material in function of the penalization parameter , is well studied and demonstrated in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF]. Borrvall and Peterson suggested starting the optimization with low penalty values and gradually modifying the penalization parameter to higher values. The function (2.15) could be written in general form:

  1 ( , ) ( ) 1
  s  f  0 1 0.01 p  0.1 p  1 p    f  s  0 1 0.01 p  0.1 p  1 p  52

Bi fluid interpolation function

The development of multi fluid interpolation function is based on interpolation functions already used in topology optimization of mechanical structure for multi material. Sigmund [START_REF] Bendsøe | Material interpolation schemes in topology optimization[END_REF] proposed a SIMP like interpolation model for three-phase interpolation, two materials and void. He stated that it is most convenient to first interpolate between the two materials and then interpolate between the resulting material and void. This formulation appears to be convenient with fluid flow and heat transfer modeling where it was assumed that fluid 1 and fluid 2 form together a single third fluid whose properties are interpolated between the two fluids, and the third resulted equivalent fluid forms a porous medium with the quantity of solid present in the cell. Following this assumption, the multi-material interpolation function introduced by Sigmund in [START_REF] Bendsøe | Material interpolation schemes in topology optimization[END_REF] appears to be convenient to the bi-fluid problem. It allows to interpolate first between the two fluids and then to interpolate between the resulting equivalent fluid and the solid. The interpolation function as expressed in [START_REF] Bendsøe | Material interpolation schemes in topology optimization[END_REF]:

  1 2 2 1 2 1 2 1 2 2 ( , ) 1
p p p E E E           (2.17)
where and are respectively the elasticity of the two materials, and the resulting elasticity of the entire cell. and are the design variables and and are their relative penalization parameters. It should be noted here that void has zero value for the elasticity, therefore no physical quantity for void phase appears in the formula. Therefore, for the case of two fluids and one solid, where the solid phase has a value for its physical quantity, the general form of the interpolation function could be expressed as follows:

    1 2 1 1 2 2 1 2 2 2 1 1 ( , ) ( , ) ( , ) 1 ( , ) 1 ( , )
f f s g p g p g p g p                    (2.18)
where , and correspond respectively to fluid 1, fluid 2 and solid physical quantities.

It is clearly shown in equation (2.18) that is used to interpolate between fluid 1 and fluid 2, and to interpolate between the resulting fluid mixture and solid. Replacing in equation (2.18) by its expression used by Borrvall and Petersson for mono fluid interpolation (equation (2.16)), equation (2.18) becomes:

1 2 2 1 1 2 1 2 1 2 2 1 1 1 2 2 2 2 1 1 1 1 1 1 ( , ) 1 1 f f s p p p p p p p p                                              (2.
        2 1 2 1 2 0, ( ) if ( ) ( ), ( ) 1, 0 if 1,1 if s f f xx x x x x x                    (2.
    12 1 2 1 2 1 2 , 1 1 2 1 2 1 2 1 2 1 2 2 min z( , ), , with z( , ) (u, , ) subject to Equations (3.4) (3.5) (3.7) 1 ( , ) 1
1 ( , ) f f f f F f P T G d N G d N                              (2.21) 54 (A) Fluid 1 constraint (B) Fluid 2 constraint

Finite volume discretization, direct problem

This paragraph presents the finite volume method applied to fluid dynamics problem as introduced by Versteeg and Malalasekera [START_REF] Versteeg | An introduction to computational fluid dynamics -The finite volume method[END_REF]. The discretisation of the physical equations is repeated to take into account the inverse permeability coefficient added to Navier-Stokes equations. Equations (2.5) and (2.7) are discretized using finite volume method over a uniform cartesian grid with control volumes at fixed size Temperature and pressure fields are evaluated on the same control volume, named scalar control volume. The velocity components are evaluated on staggered grids, where the control volume of the horizontal component of the velocity is centered on the east face of the scalar control volume and the vertical component is centered on the north face of the scalar control volume. The reason to use staggered grid for velocity is to avoid the "checker-board" problem defined by the non physical behavior when computing the pressure field if velocities are evaluated at original scalar control volumes, as discussed in [START_REF] Versteeg | An introduction to computational fluid dynamics -The finite volume method[END_REF]. Design variables and physical parameters such as density, thermal conductivity, etc are stored on the scalar grid. Equation (2.5) is projected on -axis and -axis which leads to two equations respectively for U and V, the horizontal and vertical component of velocity field, thus equations (2.5) and (2.7) are represented in generic form as follows:

    . . . .                                U V R x y x x y y (2.22)
where denotes for and and and are physical constants. They could be summarized as follows ( in bi-eta formulation):

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 ( ) u u u u u v v v v v T T P T T T U R P x V R P y T C k R q                                                                                     (2.
C C C C S S D D D D                                   (2.24)
where and represent respectively the advection and diffusive terms at cell faces and and to represent source terms, expressed as follows: .11 represents a spatial distribution of temperature, horizontal velocity component and vertical velocity component and their conduction, advection and source terms over their respective calculation grids. It should be noted that the understanding of spatial distribution of all physical quantities, design variables, physical fields and their coefficients is very important in the formulation of the discrete adjoint method, presented later in this chapter. control volume is the rectangle defined by the points A1, C1, C3 and A3 and its central node at B2, control volume is the rectangle defined by the points B1, D1, D3 and B3 and its central node at C2, and finally control volume is the rectangle defined by the points A2, C2, C4 and A4 and its central node at B3. Pressure is also evaluated at C2, similarly to the temperature .

        w e s n w e s n C U y C U y C V x C V x             (2.25)                 w e s n w e s n y y x x D D D D x x y y (2.26) Pu S x y S R x y         (2.27) P e N S E W w s n 2 x  2 x  2 x  2 y  2 y  2 y  2 y  V n V s U e U w 56
At the same point the following parameters are stored: , and . It's clearly shown that U-grid and V-grid are shifted respectively by half mesh in x-direction and half mesh in y-direction in relation to the scalar grid. and are indexes used to locate the position of each physical quantity and their terms on the corresponding calculation grid. For an optimization domain having design cells in direction and design cells in direction,. and vary as follows for every physical field grid's:

    11 11 x y U x y x y V x y x y T x y U i N j N N N N V N j N N N N T i N j N N N N                         -grid: 1 1 -grid: 1 1 -grid: 1 1 (2.28)
Where , and are the grid size of and fields. For more details on calculation of discretized terms, and considerations that have to be taken due to staggered grids, refer to [START_REF] Patankar | Numerical heat transfer and Fluid flow[END_REF].

Differencing scheme

To solve equation (2.24), transported property must be calculated at control volume faces ( at left hand side of the equation). This requires the use of an appropriate differential scheme to compute those properties in function of properties at surrounding nodes ( ):

()

face nb f   (2.29)
The choice of differential scheme becomes more critical in problems involving convection in comparison with pure diffusion problems, because transported quantity is distributed in all directions in diffusion process, whereas in convection the quantity is transported only in flow direction [START_REF] Versteeg | An introduction to computational fluid dynamics -The finite volume method[END_REF]. Available numerical schemes used are: Central differencing scheme (CDS), Upwind differencing scheme (UDS), Hybrid differentiable scheme (HDS), Quadratic upwind differencing scheme (QUICK), etc. In this work two schemes were considered: the central differencing scheme (CDS) for low Reynolds number flows and Quadratic scheme (QUICK) for higher Reynolds number flows.

Central Differencing Scheme

In central differencing scheme, the transport properties at cell faces are expressed as follows:

2 2 2 2

P W P N P S PE e w n s                    (2.30)
Replacing cell faces properties expressions above in equation (2.24) for control volume in the discretized equation is then expressed as follows (see [START_REF] Versteeg | An introduction to computational fluid dynamics -The finite volume method[END_REF] for further details):

, , , , , , , 1, 1, , 1 , 1
i j i j i j i j i j i j p i j w i j e i j s i j n i j u A A A A A S               (2.31)
Where indexes , and corresponds respectively to indexes and in equation (2.24). The coefficients of equation (2.31) are expressed as follows [START_REF] Versteeg | An introduction to computational fluid dynamics -The finite volume method[END_REF]:

, , , , , , , , , , i j
i j i j i j i j i j i j i j i j i j p w e s n e w n s p A A A A A C C C C S          , , , , , , , , , , , , 1 2 1 2 1 2 1 2
i j i j i j w w w i j i j i j e e e i j i j i j s s s

i j i j i j n n n A D C A D C A D C A D C     (2.32)

Quadratic upwind differencing scheme

Despite the advantage of CDS related to its differentiable and continuous formulation, it suffers from not taking into consideration flow direction in the computation of transport terms, and being limited to low Peclet numbers ( [START_REF] Patankar | Numerical heat transfer and Fluid flow[END_REF]. Thus for the same geometry and same control volume dimensions (which defines the number of cells) and same physical properties, increasing the velocity of the fluid requires decreasing the dimensions of the cell so the Peclet number in each cell remains less than 2. Thus increasing the velocity of the fluid requires increasing the number of cells. However the maximum number of cells in the domain is limited by the maximum number of variables the optimization algorithm could handle, in addition to the increasing in computational time when the number of cells increases. Therefore Hayase et al [START_REF] Hayase | A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures[END_REF] QUICK scheme (Quadratic Upstream Interpolation for Convective Kinetics), which stability is independent from Peclet number, is implemented to be used in topology optimization flows with high Reynolds number. Further details about QUICK scheme are presented in appendix C1.

Resolution of the equations system

After applying boundary conditions relative to every field, the problem becomes to solve the equation system made of equations ( and are the number of horizontal and vertical elements), where or . is a vector made of and is a vector made of source terms and constant corresponding to boundary conditions. is a non symmetrical pentadiagonal matrix made of coefficients (equation (2.32)) having the following form:

It should be noted that and respectively the matrices of U and V velocity component are diagonally dominant, which satisfies the boundedness property for differencing scheme. This is ensured by the presence of inverse permeability in coefficient which increases the value of . The system is solved using the Biconjugate Gradient Stabilized method [START_REF] Van Der | BI-CGSTAB: A Fast and Smoothly Converging Variant of BI-CG for the Solution of Nonsymmetric Linear Systems[END_REF].

Velocity components of equation (2.5) verify continuity equation (2.4) if the pressure field of equation (2.5) is correct. It could be seen from equations (2.4) and (2.5), that pressure and velocity variables are higher than the number of equations. Therefore, an iterative method is needed to solve the underdetermined system of equations (2.4) and (2.5) and find and fields. Thus the velocity-pressure coupling problem is solved using SIMPLER (Semi-Implicit Method for Pressure-Linked Equations Revised) [START_REF] Patankar | Numerical heat transfer and Fluid flow[END_REF], which is an improved version of SIMPLE algorithm introduced in 1972 by the same authors based on guess-and-correct iterative procedure for the calculation of pressure. A new equation system relative to pressure is introduced by SIMPLER. Finally the four systems solved in finite volume problem are: The discretized coefficients of matrix and vector are presented in details in [START_REF] Versteeg | An introduction to computational fluid dynamics -The finite volume method[END_REF].

In order to improve the convergence of the resolution procedure of the system of equations, under relaxation factor is introduced which allows equation (2.31) to take into consideration computed on previous iteration of SIMPLER iterative procedure. Equation system becomes expressed as follows [101]:

1 p UR nb nb p old nb UR UR A A B A           (2.34)
where is the under relaxation factor that varies between zero and one, the central coefficient (equation (2.31)) and corresponds to side coefficients (equation (2.32)).

Objectives functions

The total power dissipated in a porous fluid system is expressed as follows [102]:

2 2 , 1 ( , ) 2
j i i ji i j i u u u u d xx                       (2.35)
where and denotes for the directions and the velocity. Equation (2.35) evaluates power dissipation locally at each cell inside the domain. Olesen et al. [START_REF] Olesen | A high-level programming-language implementation of topology optimization applied to steady-state Navier-Stokes flow[END_REF] stated that under steady state condition and for the case where no slip condition is applied on all external solid walls and where velocity vector is perpendicular to surface at inlet and outlet, equation (2.35) could be evaluated as follows on the inlet and outlet flow sections:

2 2 11 (u, ) u u u u 22 out in d f P P d P d                       (2.36)
where denotes the outlet sections of fluid 1 and fluid 2 at the boundaries and the inlet sections of fluid 1 and fluid 2. In this work evaluation of fluid power dissipation will be evaluated by both local (equation (2.35)) and global (equation (2.36)) criterion.

Sensitivity analysis

The next step in topology optimization algorithm requires the evaluation of the objective function's derivative with respect to the design variables. It could be seen from paragraph 2.6 that depends explicitly from and hence depend implicitly from design vector ( ) in mono-eta formulation and the vector ( ) in bi-eta formulation, through the physical parameters computed by the interpolation function. Hence the gradient of the objective function cannot be computed analytically. A simple method for computing the gradient is to take a small variation in each design variable in turn and recalculate the objective function , hence the partial derivative is computed using the finite difference fomula:

( ) ( ) dF F F d        (2.37)
However the computation time can be excessive, in fact the calculation of sensitivities of variables, requires CFD calculations. A widely used technique in shape and topology optimization for computing the sensitivity of the objective function is the adjoint method, in its two forms: continuous form and discrete form. In both adjoint methods, the objective function and physical (flow and thermal) equations are combined through the use of adjoint variables. The entire combination creates the Lagrange function. The partial derivative of both physical and objective functions with respect to physical variables ( and produces the adjoint vector. Their partial derivative with respect to design variables and the adjoint vector computed produces the gradient. However the difference between continuous and discrete adjoint is the following: in continuous adjoint, all the equations are derived analytically from their initial forms (equations (2.5),(2.7) and (4.6)), and then discretized and solved. In discrete adjoint method, the same mathematical formulation is applied, but it's directly applied to the set of discretized equations (finite volume discretization method in the present case). According to Nadarajah and Jameson [START_REF] Nadarajah | A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization[END_REF], the discrete gradients have better agreement with finite difference gradients than the continuous gradients. However it was shown that this difference is small and the difference between the three methods is reduced when the mesh sized is decreased. In this work discrete adjoint approach is adopted for the calculation of the gradient, and is presented in details below.

Discrete adjoint approach

Consider an objective function subject to , the residual of the resolution of system of equations in finite volume method. The Lagrange function is introduced and defined as follows [START_REF] Giles | An Introduction to the Adjoint Approach to Design[END_REF]:

( , ) ( ) . ( , ) t L X f X J X     (2.38)
where is the adjoint vector, with the same dimension as the vector. When satisfies the Lagrange function is equal to the objective function The total derivative of equation (2.38) with respect to is:

t dL df dJ d d d      (2.39)
Equation (2.39) could be written as follows in function of partial derivatives of the objective and residual functions:

tt dL J f f J X d X X                       (2.40)
The critical term to evaluate in equation (2.40) is . Discrete adjoint approach suggests that the adjoint vector verifies the equation multiplied by , hence it doesn't need to be evaluated:

0 tt t f J J f X X X X                            (2.41)
Finally after computing the adjoint vector using equation (2.41), known as adjoint equation, the total derivative of the objective function is established by the following equation:

t df J f d        (2.42)

Calculation of  df d

The direct variable vector , design variable and adjoint vector are expressed as follows:

 

1 2 = = = U V mono eta bi eta P U XV P                          (2.43)
The residual function , which should verify , is the discretized system of equations of finite volume method: . ( , ) .

.

U U U V V V P P P J A U B J X J A V B J A P B                 (2.44)
The evaluation of aims at solving the following systems of equations:

tt t Jf XX df J f d                               (2.45) With matrix J X   expressed as follows: 63 UUU VVV PPP JJJ U V P JJJ J X U V P JJJ U V P                         (2.46)
And the vector J    is expressed as follows for mono eta and bi -eta fomulations:

12 12 12 U UU V V V PP P bi eta mono eta J JJ J J J JJ JJ J                                                 (2.47)
Finally the partial derivatives of the objective functions are the following:

1 2 mono eta bi eta f f U f f f f f XV f f P                                                 (2.48)
. is a vector which elements are evaluated by deriving equation (2.36) with respect to and over their calculation grid. The matrix has the same number of rows in monoeta and bi-eta formulations, but the number of columns doubles in case of bi-eta formulation because each cell holds two design variables. In next paragraph details about the calculation of , and are presented. The assembly of , and lies on the derivatives of with respect to every term of a vector where or and or . The derivative of with respect to is expressed as follows:

X X X X J A B X XA                (2.49)
It should be noted that if and 0 otherwise. The critical term of equation (2.49) is the evaluation of matrix , which requires first to evaluate the derivative of the coefficients and of field for and with respect to , ( and ) )then to evaluate the derivative of interfacial and center coefficients.. The derivatives of coefficients relative to central differencing scheme are:

, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , 1 2 1
A S A A A A C C C C A D C A D C A                                                              ,, , , , , , , , , , , , ,, 1 2 1
2 a b a b ss i j i j i j a b a b a b n n n i j i j i j ab ab u i j i j DC A D C S B                                             (2.50)
The same method is applied for QUICK scheme, but using the derivative of its corresponding interfacial and center coefficients, presented in Appendix C1. The matrix is constructed in a similar manner to the matrix in finite volume method. Thus for every element of the following set of equations are computed:

0 0 , , , , , , , ,
.

X X X i j i j i j X X i j i j X X X i j i j i j r r A A B X A B X A A B X                                                                                               (2.51)
The result of (2.51) is a vector which dimension is equal to the number of elements of grid.

Finally the matrix X J    is constructed as follows:

(2.52)

With and For or and for or P , the sub matrices are evaluated and then the big matrix is assembled as seen in equation (2.46). Simlarly, for or and for the sub matrices that computes derivative of discritized finite volume equations with respect to design variables are evaluated and the matrix is assembled as seen in equation (2.47). For a domain having design variables, is a square matrix with , and is a rectangular matrix for mono eta formulation and for bi fluid formulation with . Finally it should be noted that a special attention should be given to the spatial distribution of all finite volume coefficients, physical and design parameters.

Regularization techniques

In topology optimization the following numerical problems could be encountered: checkerboard that happens when solid and void appear in adjacent meshes (or solid and fluid), mesh dependency when different results are obtained when using different mesh sizes for the same problem and the appearance of intermediate state or porous cells in the final design [START_REF] Sigmund | Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima[END_REF]. Regularization techniques are used in order to avoid the checkerboard problem.

In their work [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF], Borrvall and Petersson concluded that there is no need for regularization techniques in fluid mechanics topology optimization problems for solving checkerboard. On the contrary Wiker et al. [START_REF] Wiker | Topology optimization of regions of Darcy and Stokes flow[END_REF] proved that regularization techniques are needed in such problems. However, filtering can still be applied in order to ensure that no checkerboard will appear in the final solution. Filtering methods for topology optimization are divided into two categories, density based and sensitivity based. Density based filters [START_REF] Bruns | Topology optimization of non-linear elastic structures and compliant mechanisms[END_REF], [START_REF] Bourdin | Filters in topology optimization[END_REF] are directly applied to the

X N J                      X N  ,, . XX i j i j AB X      0 X A if X if X        
design variables, while sensitivity filters [START_REF] Sigmund | Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima[END_REF], [START_REF] Sigmund | Design of Material Structures Using Topology Optimization Phd Thesis Department of solid mechanics technical university of denmark[END_REF] are applied to the computed sensitivities. Density filters consist of replacing the design variable in each element of the domain by a weighted average of the elements in a mesh-independent neighborhood of this element. It is applied to the design variable in each iteration after the optimal solution is recalculated by the optimizer and before CFD solution is calculated. Sensitivity filter on the other hand, is the redefinition of the sensitivity at an element as a weighted-average of the sensitivities in a meshindependent neighborhood of this element. The neighborhood of an element is defined in general as the elements whose centers are in a pre-defined distance from the center of this element. In this work both sensitivity and design filters were used. Density filter of Bruns et al. [START_REF] Bruns | Topology optimization of non-linear elastic structures and compliant mechanisms[END_REF] was considered, which is defined as follows:

1 ˆN i ij i j w      (2.53)
where is a weighting factor, is the new density after applying the filter at i-th node and the initial density value at the same node. It exist in literature a wide range of weighting factors applied to filters, as linear weighting, Gaussian weighting, etc. In this work the following weight factor is used [START_REF] Svanberg | Density filters for topology optimization based on the geometric harmonic means[END_REF]:

ij ij il l rd w rd     (2.54)
where is the filter radius defined by the user and is the distance between the central nodes of and cells, where is a cell that exist within a circle of radius centered at the central node of cell .

is equal to 0 for elements outside the circle . is the distance between and elements, where denotes for all elements within the circle . It should be noted that the choice of the filter radius is very delicate and has a major impact on the convergence to the optimal solution. Usually the filter radius varies dynamically throughout the optimization and it should have a low value at the final iterations to allow the algorithm to reach black and white solution. Using the same weighting factor, the following filter is applied on sensitivity field [START_REF] Svanberg | Density filters for topology optimization based on the geometric harmonic means[END_REF]:

1 ˆ1 N ij j j i i j df df w d w d      (2.55)

Optimizer

The Method of Moving Asymptotes (MMA), adopted in this work, was developed by Svanberg in 1987. The method is based on a special type of convex approximation [110] [START_REF] Svanberg | The method of moving asymptotes-a new method for structural optimization[END_REF]. It can fit single and multiple constraints, and proved to be efficient in topology optimization problems. Svanberg proposed a better convergent version of MMA in 2002, the globally convergent version of the method of moving asymptotes (GCMMA) [START_REF] Krister Svanberg | A class of globally convergent optimization methods based on conservative complex separable approximations[END_REF]. However GCMMA needs relatively higher calculation time. MMA has become a standard optimization tool for topology optimization problems, especially with density based methods [START_REF] Deaton | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF].

MMA solves the optimization problem as follows: suppose having the following optimization problem:

0 ,min ,max minimize ( ) subject to ( ) 0 1,..., 1,..., k k ii k k k j j j f f b i m jn            (2.56)
where is the number of constraints, the number of variables and the index of the iteration in the optimization loop.

At each iteration a subproblem is created, in which the functions (objective function and constraint functions) are replaced by an approximated convex function [START_REF] Svanberg | The method of moving asymptotes-a new method for structural optimization[END_REF]. The approximated functions are computed on basis of the gradient information of the current iteration point and the moving asymptotes of the current iteration. The moving asymptotes are the lower and upper asymptotes in function of current iteration point , the previous iteration point and the asymptotes of the previous iteration, and a dimensionless parameter to limit the search area of optimal point between the asymptotes. Many authors proposed new methods to update the approximation functions [START_REF] Deaton | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF] [114], and to improve the convergence of the original method proposed by Svanberg [START_REF] Svanberg | The method of moving asymptotes-a new method for structural optimization[END_REF]. The subproblem generated by MMA is then solved using "primal-dual-interior-point method" to find the optimal point that will be used in next iteration.

Results

Mono-eta formulation

The configuration presented in Figure 2.12. is a square domain of 36 side length having an inlet and outlet section for each fluid respectively at west and east sides The maximum porosity for each fluid is is , which is equivalent to the volume of fluid required by a straight line between the inlet and outlet of each fluid. Both fluids enter and leave the domain with same velocity profile:

. The inverse permeability coefficients are: and . The remaining of fluid physical properties are considered constant in equation (2.5) and . This means that only is interpolated using equation (2.12) The size of the design grid is . Four different simulations were conducted in the present case study, which initial configuration is depicted in Figure 2.13. The four simulations differ by the initial guess of field, in order to analyze the behavior of mono-eta interpolation function regarding its stationary point, as discussed in the previous paragraph. The first column in Figure 2.13 shows the homogeneous initial domains whereas the second column shows their corresponding topology optimization results ).The color map used to represent the values of design variables field in the domain varies from blue to grey passing by light blue (0 to 0.5) and from grey to red passing by light red (0.5 to 1). It could be clearly seen from topology optimization results in Figure 2.13 that approximately all cells converged to one of the three phases: fluid 1, solid or fluid 2. As discussed before, the convergence to pure solid, pure fluid 1 or pure fluid 2 cells( equal to 0, 0.5 or 1) has some requirements on the shape of the interpolation function. It has been found that is required to ensure the total convergence to the discrete values that defines the three phases. However despite reaching 0-0.5-1 discrete solution with , it was not sufficient to reach the global optimum and further increasing of was required. Finally the optimal solution seen in Figure 2.13.D was reached with . For this value of , inverse permeability coefficient at and is much higher than , which will lead to unrealistic and underestimated values in velocity fields. Hence the values of the objective function in Table 2.1 are much higher than the real values of fluid power dissipation caused by fluid motion. Regarding the final optimal shape reached in function of initial value of , results showed:

 For , the optimization algorithm was not able to add nor fluid 1 nor fluid 2 in the domain and all cells remained pure solid, (Figure 2.13.A). The reason is the null derivative of the objective function at this point.

 For

, the algorithm begins to add only fluid 2, represented by red color, between inlet and outlet sections of fluids. However the two pipes formed by fluid 2, don't fill the complete section between the inlet and outlet boundaries of each pipe, due to the porosity constraint of fluid 2 that is equal only to the volume of one complete pipe. Fluid 1 quantity represented by blue color, which was added at later stages of optimization process couldn't be added in the optimum place to complete fluid pipes which show the difficulty to add fluid 1 when the initial value of design variable is greater than 0.5, as seen in Figure 2.14. It should be noted that same structure is reached when the initial value of is higher than 0.51.

 For

the same behaviour as for is observed but with convergence towards fluid 1 instead of fluid 2.

 For the case where in the upper half of the domain and in the lower half of the domain, the algorithm was able to reach global optimal solution and to optimally distribute fluid 2 in the region where and fluid 1 in the region where . It could be clearly seen from Figure 2.13.D that the allowable porosity of each fluid occupied the straight line between inlet and outlet section, which is the optimum solution for this kind of configuration. This The initial configuration for this case is the same presented in Figure 2.12, but with side length , and a design grid size . The maximum porosities are . The objective function minimized is the global fluid power dissipation (

). Figure 2.15 shows the final result of field, where the black color represents solid cells ( ) and white color represents fluid cells (

). It should be noted here that white color cells are made of 50 % fluid 1 and 50 % fluid 2, which means a mix of both fluids. The design variables converged to their upper and lower limits, known as black-and-white solution, whereas remains equal to 0.5, its initial value, in all design cells. Density and dynamic viscosity are constants in equation (2.5) ( and ), whereas is computed using equation (2.19), with and .

Topology optimization result of bi-fluid problem using bi-eta interpolation function depends on the convergence of two design variables in each cell:

for the separation of solid and equivalent fluid phases, and for the separation of fluid 1 and fluid 2 phases in the fluid portion of the cell. It is obvious here that converged to a discrete solution, whereas remained equal to its initial value in all design cells. This behavior was predictable because the term is constant, hence equation (2.19) is similar to (2.15), which means the gradients of the objective and constraint functions with respect to are null ( and ) Thus separation of fluid 1 and fluid 2 phases has no effect on the minimization of the objective function nor on constraint functions.

Regarding the convergence of , it could be seen that the quantity of each fluid occupy the straight line between the inlet and outlet sections, which is equal to the maximum volume of fluid allowed by the porosity constraints. This means that the final solution reached by the algorithm is the real optimal solution which is known a priori for this configuration. Figure 2.16 shows the minimization of the power dissipation function. At first the objective function decreases quickly while the algorithm tends to reallocate the fluid from the entire domain to form a continuous fluid flow between the inlet and outlet boundaries. After iteration 24, an approximated architecture for the channel is visible, and the algorithm changes the values of design variable at the solid fluid interfaces to get the optimal shape of Figure 2.15. Figure 2.17 shows that ( ) continue decreasing after iteration 24, to reach the optimal solution after 61 iterations approximately. It should be noted also that sensitivity filter was used at the beginning of the optimization to help the algorithm avoiding checkerboard problem. This filter was omitted later to allow the convergence to a discrete 0-1 solution for field.

Single pipe example: non-symmetric properties for fluids

In the second example, an imbalance in fluid 1 and fluid 2 properties and maximum porosities is created to study the convergence of design field. In the previous example inverse permeability coefficients of both fluids and where equal to 0, and the maximum allowable porosities of each fluid in the domain and were also equal, remained equal to 0.5 in all cells, hence all fluid cells were made of a mixture of 50 % fluid 1 and 50 % fluid 2.

The configuration presented in Figure 2.18 consists of a single inlet and single outlet for fluid flow located respectively at west and east sides of the domain. The remaining boundaries are adiabatic non slip walls. The Reynolds number of the fluid based on the average inlet velocity and its inlet dimension is . The dimension of design variable grid is . The color map is defined as follows: when the values go from 0 to 1, representative color goes from black, to dark red, then red, then yellow, then light yellow and finally to white. These three examples shows the possibility of the convergence of and fields simultaneously, hence to optimally distribute two fluids and one solid phase in a fixed domain using the bi-eta penalized interpolation function developed in paragraph 2.4.2.

Case studies

As seen in last paragraph, bi-fluid topology optimization problem is solved similarly to mono-fluid problem in bi-eta formulation when the two fluids have similar physical properties and maximum porosity limit. In this paragraph, identical properties and porosities limits for fluid 1 and fluid 2 will be considered, which allow comparing the results of the method with mono fluid literature results. Evaluation of fluid power dissipation by global and local objective functions (respectively equations (2.35(2.36) will be considered. CDS and QUICK schemes and their corresponding discrete adjoint method will be also considered. It should be noted that in previous results global objective function and CDS scheme were used.

Diffuser

.21: Initial configuration of diffuser case

The diffuser case depicted in Figure 2.21 was studied in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] and [START_REF] Guest | Topology optimization of creeping fluid flows using a Darcy-Stokes finite element[END_REF] for low Reynolds numbers. The fluid enters the square shaped domain from the whole west edge, and leaves it from a section located at the center of the east edge and covering one third of its length. To ensure mass continuity through the domain, the velocity at outlet section is three times bigger than the inlet velocity. The density of both fluids is and the dynamic viscosity . The maximum porosities of fluid 1 and fluid 2 are respectively and . have an initial value in the whole domain equal to 0.5, and because of the equality of fluids properties and maximum porosities, remains constant throughout the optimization. Hence the problem is similar to a mono fluid optimization problem with an equivalent fluid having the same density and dynamic viscosity as fluid 1 and fluid 2, and having a maximum porosity limit . The inverse permeability coefficient is equal to 0 for the equivalent Global and local fluid power dissipation are considered in this example. For global power dissipation function, the velocity components and belongs to the boundaries outside the design domain, which means that the partial derivative of the global power dissipation function with respect to velocity is null. Pressure components are shifted by half mesh from velocity components, hence they are evaluated at meshes inside the design domain facing fluid inlet and outlet boundaries. Figure 2.22 shows an example of the spatial distribution of pressure and velocity components involved in the calculation of global power dissipation. and are indexes to locate the meshes at inlet and outlet velocity boundary conditions. It should be noted also that the global power dissipation function doesn't depend directly from design parameter, hence its partial derivative with respect to is null. Equation (2.48) for global power dissipation is then expressed as follows (

):

0 0 0 0 dd bi eta d ff X f P                   (2.58)
4 i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j U U U U U xx U U U U U yy V V V V V xx V V V V V yy                                           (2.59) 1, , 2, , 1 , 1 1, 1 1, 1 1, 1, 1, 1 1, 1 , 1 , 1 , , 2 4 4 4
The presence of inverse permeability coefficient in the expression of local fluid power dissipation function makes it directly dependent from design variables; hence the equation(2.48) is expressed as follows (

):

1 2 0 bi eta U XV                                  (2.60)
Optimization results using both objective functions are presented in Figure 2.24. The simulations were conducted using central differencing scheme. The only difference that could be seen between the two structures is that in case of local evaluation of objective function, the fluid

1, 1 i j U   , 1 ij U  1, i j U  , ij U 1, 1 i j U   , 1 ij U  , ij  , ij  x  1, i j U  2, i j U  , ij V , 1 ij V  , 1 ij V  , 2 ij V  1, i j V  1, i j V  1, 1 i j V   1, 1 i j V   y 
channel's upper and lower limits at west side are located at the north west and south west corners of the domain, whereas for global objective function these limits are a little bit shifted toward the center of the west edge. However this could be considered as a minor difference which has not a significance impact on fluid power dissipation. Hence both objective functions had led to the same structures which are similar to the result of Borrvall and Petersson in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF]. Finally it should be noted that the number of iterations to reach the final optimal solution for both simulations were similar too (72 iterations for the global objective function and 78 iterations for the local objective function). The second case is the bend pipe pictured in Figure 2.25 which was considered in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] for low Reynolds numbers and in [START_REF] Deng | Topology optimization of steady and unsteady incompressible Navier-Stokes flows driven by body forces[END_REF], [START_REF] Gersborg-Hansen | Topology optimization of channel flow problems[END_REF] and [START_REF] Sigmund | Topology optimization for multiphysics problems: A future FEMLAB application?[END_REF] for low to moderate Reynolds number. The domain is square shaped with side length . The fluid enters the domain from the west side and leaves it from the south. The same equivalent fluid properties of diffuser case are considered in this example. The maximum equivalent fluid porosity is . In this example QUICK scheme was used in finite volume method with its corresponding discrete adjoint in sensitivity analysis, which allows to consider higher Reynolds number than in the case of using CDS scheme. Fluid power dissipation was evaluated using the global criteria. We can see that for low Reynolds the pipe is straight, whereas for the pipe has an arc shape. The same shape variation of the bend pipe with respect to the increase of Reynolds number were observed in [START_REF] Deng | Topology optimization of steady and unsteady incompressible Navier-Stokes flows driven by body forces[END_REF], [START_REF] Gersborg-Hansen | Topology optimization of channel flow problems[END_REF] and [START_REF] Sigmund | Topology optimization for multiphysics problems: A future FEMLAB application?[END_REF]. Sigmund [START_REF] Sigmund | Topology optimization for multiphysics problems: A future FEMLAB application?[END_REF]. stated that the bend has sharp corners for low Reynolds number, but the corners become rounder when the inertia term increases, which are also observed in this study. To investigate the topology optimization results, fluid power dissipation was computed for and for the structures seen in Figure 2.26.A, B and C. Table 2.2 shows an analogy between the results of topology optimization in Figure 2.26 and the results of the investigation conducted.

Conclusion

In this chapter the method of bi-fluid topology optimization was introduced, with detailed explanation of each part of the algorithm. The method is based on coupling the Finite Volume Method for the direct solver, the discrete adjoint method for sensitivity analysis and the method of moving asymptotes for the non linear optimization problem. In order to introduce threephase (two fluids and one solid) to the fluid flow problems, two interpolation functions for density approach were developed, implemented an tested : mono-eta interpolation function (using a single design variable in each cell), and bi-eta interpolation function (using two design variables and in each cell).

Mono-eta formulation has the advantage that it uses fewer variables than the bi-eta formulation but suffers from major drawbacks that limit its use. First the convergence shows high dependency to initial guess of the design field

, where the convergence of the algorithm to the global optimal solution was possible only when the initial design field is set up with convenient values. Second the convergence to global optimal solution needs high values of penalization parameter , in this case the permeability coefficients of the two fluids are not equal to 0, which affect the modelling of fluid phases in Navier Stokes equation. Despite these inconvenient mono-eta interpolations function stills an interesting formulation, especially when higher number of phases has to be handled.

In Bi-eta interpolation function each design cell holds two design variables, to interpolate between solid and fluids phases, and to interpolate between fluid 1 and fluid 2 phases. Results showed the capability of convergence of both design fields simultaneously to optimal solution without being affected by the initial guess, hence this method optimally distribute one solid phase and two fluids phases only in function of minimizing the objective function subject to constraint functions and physical equations of the problem.

For designing bi-fluid engineering devices using bi-fluid topology optimization an important feature should be considered: separation of fluid phases under any set of parameter values and forcing each fluid to take its own path throughout the domain without crossing or mixing with the channel of the second fluid. Fluid 1 channel should consist only of 100 % fluid 1 cells, and similarly fluid 2 channels should be made only of 100 % fluid 2. This will be the topic of the next chapter, in which many methods will be implemented to force the separation of fluids subdomains, avoid their intersection and preserve a minimal solid thickness (at least one solid cell) between them. 3.1 montre que la fonction de continuité n'aura pas des valeurs positives sauf dans le cas où deux fluides différents sont présents dans deux mailles voisines. La fonction multi-objectif, couplant la fonction de perte de charge et la nouvelle fonction de continuité, est exprimée dans (3.6).

Fonction objectif de continuité a) Formulation et implémentation

b) Résultats

Le problème d'optimisation consiste maintenant à minimiser la fonction de perte de charge en assurant que la fonction de continuité converge vers 0 dans la solution finale Figure 3

Modification du coefficient inverse de perméabilité.

Dans le paragraphe précédent, les sous-domaines fluide 1 et fluide 2 sont séparés en utilisant une fonction objectif qui pénalise la présence de deux fluides différents dans des mailles voisines. Cette fonction objectif a été ajoutée à la fonction multi-objectif unique minimisée par l'optimiseur. Dans ce paragraphe, une autre méthode est formulée pour forcer la séparation des sous-domaines fluides, en agissant indirectement sur la fonction de dissipation de puissance du fluide par l'intermédiaire du coefficient de perméabilité inverse dans l'équation de Navier-Stokes. Le coefficient de perméabilité initial ralentit la vitesse d'écoulement proportionnellement à la quantité de solide dans la maille. Cependant, un autre terme proportionnel à la vitesse du fluide est ajouté à l'équation , et agit comme une force interne ou un terme d'absorption au champ de vitesse lorsque deux fluides différents sont présents dans la même cellule ou des cellules adjacentes. Ceci pourrait être vu comme si le fluide 1 agissait comme un mur pour le fluide 2 et vice versa, ce qui augmenterait la dissipation de puissance totale du fluide. La nouvelle expression du coefficient inverse perméabilité est présentée dans les équations (3.8) à (3.13).

Les cas d'études des doubles tuyaux et des doubles diffuseurs sont reconsidérés dans la validation de la 2 ème méthode de séparation. Dans le cas des doubles tuyaux, la géométrie des canaux et le positionnement des deux fluides sont similaires aux résultats de l'optimisation utilisant la première méthode de séparation. Toutefois, pour le cas des doubles tuyaux, la séparation des canaux des différents fluides est assurée, mais la forme des canaux n'est pas de la même efficacité que celle obtenue par la première méthode de séparation des fluides.

Fonction contrainte a) Formulation et implémentation

La fonction de continuité (3.1), est considérée ici comme une fonction contrainte au lieu d'être intégrée dans la fonction multi-objectif comme dans la première méthode de séparation. Le nouveau problème d'optimisation est défini dans (3.15). est un petit nombre, par exemple , pour forcer la fonction contrainte de continuité à converger vers 0 et ainsi assurer la séparation des fluides.

b) Résultats

La configuration initiale est représentée dans la Figure 3.13, et les résultats d'optimisation topologique dans la Figure 3.14. Les résultats montrent que a convergé vers 0 et 1 alors que a convergé vers 0.95 et 0.05, c.à.d. que la séparation totale des deux fluides n'est pas assurée. De plus, la méthode a besoin d'un grand nombre d'itérations en comparaison avec les deux premières méthodes de séparation. Enfin, le processus d'optimisation montre que la séparation des fluides ne commence qu'après que l'architecture des canaux soit connue, ce qui signifie dans la configuration des doubles diffuseurs que la méthode de séparation par fonction contrainte ne pourra pas empêcher la fusion des deux tuyaux au milieu du domaine. Bi-fluid topology optimization method presented in chapter 2 still lacks to one of its most important features: its ability to separate the two fluids; hence it cannot be considered bi-fluid topology optimization method yet. Paragraph 2.10 shows that the problem is still solved similarly to mono-fluid topology optimization, where the fluid present in the domain is a mixture of fluid 1 and fluid 2 (equivalent fluid). The objective of this chapter is to find the most suitable method to:

Cas d'étude

 First to prevent fluid 1 and fluid 2 from mixing in the same cell when the algorithm converges to the final solution.  Second to ensure the continuity of each fluid, which means that each fluid is transported in a single or multiple pipes that connect inlet and outlet sections of the corresponding fluid, which are predefined by the user on domain's boundaries.  Third to ensure that the channels of two different fluids should not mix or cross, and the algorithm should be able to maintain a minimal solid thickness between the channels.

For this purpose, three different methods were suggested in this chapter; hence fluid separation can be achieved:

1. By introducing a new objective function that penalizes the existence of two fluids in the same cell, or neighboring cells. 2. By taking the penalization effect as a constrain function. [START_REF]The Independent[END_REF]. By modifying the inverse permeability coefficient (paragraph 2.2.1) where fluid 1 will be an obstacle for fluid 2 and vice versa, which will affect directly the total fluid power dissipation objective function .

In order to select the most suitable technique, the three methods were tested and compared for various case studies.

Continuity objective function

Implementation

The optimization method starts with an initial domain in which all cells have the same material composition, made of a mixing of the two fluids and the solid. Progressively throughout the optimization process the method should completely separates the two fluids and the solid phases. The separation of both fluids phases from solid phase is acheived thanks to the presence of the inverse permeability coefficient in Navier-Stokes equation and to its penalized interpolation function, while there are no physical or mathematical reasons to separates fluid 1 and fluid 2 phases or to prevent them from existing in the same cell or in two neighboring cells. Therefore an objective function is formulated and added to the multi-objective optimization problem in order to mathematically prevent the mixing. The objective function, called continuity function , for a domain having horizontal meshes and vertical meshes is expressed as follows:

, ,

, xy NN ij ij c xy f f NN    (3.1)
Figure 3.1: (i,j) Mesh and its neighbor meshes where is computed in each cell, in function of fluid 1 and fluid 2 quantities in mesh and its neighboring meshes, sketched in Figure 3.1 as follows:

, , 1, 1, , 1 , 1 , 1, 1,
, 1 , 1 1 2 2 2 2 2 1 1 1 1
11 44

i j i j i j i j i j i j i j i j i j i j i j f f f f f f f f f f f Q Q Q Q Q Q Q Q Q Q                         (3.2)
where and are respectively the quantity of fluid 1 and fluid2, expressed as follows in mesh:

i,j+1 i,j i,j-1 i+1,j i-1,j   , , , 1 1 2 1 i j i j i j f Q   (3.3) , , , 2 1 2 i j i j i j f Q   (3.4)
Table 3.1 shows the values of for mesh in function of all possible combinations of phases in a particular cell ( or and or ) and its neighbors, as seen in Figure 3.1.

Mesh (i,j) Neighboring meshes

A 1 A 2 B 1 B 2 f i,j Fluid 1 Fluid 1 1 0 0 1 0 Fluid 1 Fluid 2 1 1 0 0 1 Fluid 1 Solid 1 0 0 0 0 Fluid 2 Fluid 1 0 0 1 1 1 Fluid 2 Fluid 2 0 1 1 0 0 Fluid 2 Solid 0 0 1 0 0 Solid Fluid 1 0 0 0 1 0 Solid Fluid 2 0 1 0 0 0 Solid Solid 0 0 0 0 0
i j i j i j i j i j i j i j i j i j i j i j f f f f f f f f f f f Q Q Q Q Q Q Q Q Q Q AB AB                           (3.5)
It could be clearly seen from Table 3.1 that the continuity objective function takes values greater than 0 when the cell and is neighbors are made from different fluids; hence the objective function will only prevent two adjacent cells from containing two different fluids. It should be noted also that equation (3.5) prevents also that two fluids exist simultaneously in the same cell when the algorithm converges to the final optimum.solution Finally continuity function is minimized simultaneously with the minimization of fluid power dissipation function using a multi-objective optimization. Hence the multi-objective function is expressed as follows:

d c c F f w f  (3.6)
where is a constant weight to control the preference of continuity function in the mutliobjective function. Continuity function depends directly from design parameter, hence its derivative with respect to is computed analytically. The gradient of the multi objective function with respect to design variable is then expressed as follows:

dc c df df dF w d d d     (3.7)

Results

Double pipe

The double pipe example presented in Figure 3.2.A, already considered in chapter two, is reconsidered here, with the implementation of continuity objective function in the multiobjective optimization process. The final result is shown in Figure 3. It should be noted that attention must be given here to careful choice of the values of all numerical parameters and filters applied, to ensure the convergence of the problem. We mention the most critical ones:  Numerical simulations showed that the best suitable value of the ratio is around 0.1 (at the exception of first 20% of the iterations where drops quickly due to the reallocation of the fluid). However this could slightly vary for a different configuration. It should be noted that choosing a small value of continuity function weight will not lead to fluid separation, or it needs an extremely high number of iterations for separating fluids subdomains. In contrary choosing a big value leads to a failure in the resolution of the main physical problem of minimization of the total pressure drop and creation of fluid channels.  A density filter must be applied on field, to prevent the convergence to a local minimum, which leads to fluid 1 and fluid 2 cell in the same fluid channel. The filter radius depends on the mesh size, and is progressively reduced throughout the optimization, to be completely omitted at final iteration to permit the total convergence of values to 0 and 1. Figure 3.3 and Figure 3.4 show respectively the variation of fluid power dissipation function and the total multi objective function that includes the sum of and . It could be clearly seen that the shape of multi objective function curve is similar to power dissipation objective function curve, because of the weight that has to be relatively small, hence power dissipation objective function is the dominant term in the multi objective function. Finally it should be noted that 105 iterations were required to reach the global optimal solution, whereas in the case of minimization of fluid power dissipation alone 61 iterations were sufficient, as shown in Figure 2.16. It should be mentioned that at iteration 40 the shape of channels is approximately known and the remaining of iterations aim to adjust fluid-solid interface and to fully converge to a discrete 0-1 solution. Previous example showed that mono-fluid and bi -fluid topology optimization (with and without continuity objective function) leads to the same device architecture, the only difference is the nature of fluid filled in the channels. However this is not the case for all configurations, where for certain domain geometries a contradiction between fluid power dissipation and continuity objective functions may occur, which will lead to different fluid channels architecture. It could be clearly seen that for the first case in which continuity objective function is not considered, the two fluids form a single equivalent fluid and takes a single path inside the domain by joining the pipes in the middle and disjoining them near outlet region. However, for the second case where continuity objective function is considered, each fluid takes a separate path between its corresponding inlet and outlet sections, thus two different geometries are generated. For the first case total fluid power dissipation is , whereas for the second case , which justify the fusion of both fluid channels in the first case. However with the insertion of continuity function, the minimum of the global multi-objective function corresponds to the structure that avoids pipes fusion. This example showed that using continuity objective function as a penalty function will ensure the total separation of fluids and omit all intermediate values of . The algorithm will search for the distribution of fluid channels to have the minimal pressure drop but with respecting the condition of fluid separation

Modification of the inverse permeability coefficient 3.3.1. Implementation

Inverse permeability is added to Navier-Stokes equation, to act like a friction force between the fluid and the solid in cells having intermediate values of . The shape of the interpolation function, ensures the convergence of to 0-1 solution (or near 0 and near 1 solution). Using this form of inverse permeability coefficient will lead to an optimal convergence of field, without having any effect on field. In fact, the inverse permeability coefficient of fluid 1 and fluid 2 are both equal to zero , the presence of fluid 1, fluid 2 or any mixture of them will have no effect on velocity and pressure fields, thus no effect on power dissipation objective function. Therefore the gradient of power dissipation objective function with respect to is null.

 Fluid 1 Fluid 2 u 1 1 f C  Fluid 2 Fluid 1 u X 2 2 f C  d f 95
In the previous paragraph fluids subdomains are separated using an objective function that penalizes the presence of two different fluids in the same cell or in two neighbor cells. This objective function was added to the single multi-objective function minimized by the optimizer. In this paragraph another method is formulated to force the separation of fluids subdomains, by acting indirectly on the fluid power dissipation function throughout the inverse permeability coefficient in Navier-Stokes equation. Similarly to the initial expression of that slows down the velocity proportionally to the quantity of solid in each design cell, another term proportional also to fluid velocity, is added to Navier-Stokes (equation (2.5)), and acts like an internal force or absorption term to velocity field when two different fluids are present in the same cell or adjacent cells. This could be seen as if fluid 1 acting as a wall for fluid 2 and vice versa, which will increases the total power dissipation.

The new expression for the inverse permeability coefficient is expressed as follows:

, , , , , , 1 1 2 2
i j i j i j i j i j i j f f X C C       (3.8)
where the coefficients are defined as follows:

  ,, 1 , 1 1
i j i j s f s ij p p            (3.9)
     ,, 1 1 2 1 2 , 2 1 1 1
i j i j f f f ij p p               (3.10)   ,, 2 1 2 1 2 , 2
1

i j i j f f f ij p p              (3.11) , , 1, 1, , 1 , 1 1 1 2 2 2 2
i j i j i j i j i j i j f f f f f f C Q Q Q Q Q           (3.12) , , 1, 1, , 1 , 1 2 2 1 1 1 1
i j i j i j i j i j i j f f f f f f C Q Q Q Q Q           (3.13)
and are respectively the quantity of fluid 1 and fluid 2 inside the cell, defined in equations (3.3) and (3.4) in function of and . The penalization parameters and are used to adjust the convexity of the functions as pictured in Figure 3.9. Navier-Stokes equation (2.5) becomes: The lower limits and are usually taken equal to 0, whereas the choice of the values of the upper limits and controls the relative importance of each term in the expression of .

  2 1 2 1 2 1 2 ( , ) u. u ( , )u ( , ) u XP               (3.

If

, then fluid solid interface will leads to bigger increasing in power dissipation function than fluid 1-fluid 2 interfaces. In contrary if then two neighbor cells having different fluids will have a major negative impact on fluid power dissipation objective function than solid-

1   s  f  0 1 0.01 p  0.1 p  1 p  2  1  2 f  1 f  0 1 0.01 p   0.1 p   1 p   2  2  2 f  1 f  0 1 0.01 p   0.1 p   1 p  
97 fluid 1 or solid-fluid 2 interfaces or solid and any fluid mixture interfaces. Finally a special attention must be taken in the calculation of the gradient of the objective function using the discrete adjoint method when using the new inverse permeability coefficient. The reason is that the new coefficient depends on design variables of the neighbors cells and not only the cell itself as in the initial coefficient. and . It could be seen that , which means that a preference was given to channels architecture optimization over the optimization of the nature of fluid flowing in each channel, otherwise no fluid channels will be created. However the problem converged progressively to the optimal solution by creating two straight pipes and totally separating the cold and hot streams, similarly to example 3.2.2.1. The final solution is shown in Figure 3.10 A. In comparison with the method of fluid separation using continuity objective function (paragraph 3.2), this method was able to converge to the same global optimal solution. ). It could be clearly seen that the fluid power dissipation function increases dramatically by 37517 %. This explains why the algorithm will begin to separate the fluids from the first stages of the optimization process by converging the values of to 1 near the inlet and outlet sections of fluid 2 and to 0 near the inlets and outlets sections of fluid 1, and later along the entire pipes as long as they are created progressively by the convergence of field. However small convergence steps in MMA method and a density filter on field are applied to avoid the convergence to local optimal solution similar to Figure 3.10.C. ) used in double pipe example in paragraph 2.10 and the actual inverse permeability coefficient used in this example to separate fluids subdomains. It should be noted that in both cases, the domain is initialized with the same fields of and . For the initial configuration, in case of using the inverse permeability coefficient of equation (3.8), whereas when the physical problem is solved with the initial inverse permeability coefficient presented in chapter 2. The big difference in the values of both cases is caused by the presence of a fluid mixture in all cells in the initial configuration. Hence, the flow is disturbed in case of using the coefficient of equation (3.8), because fluid 1 and fluid 2 are present simultaneously in all design cells. However both curves converge to the same minimal value of by reaching the same shape for fluid channels. It should be noted that the effect of fluid 1-fluid 2 interfaces on velocity and pressure fields using inverse permeability coefficient reformulated is totally omitted as the two subdomains are completely separated at the end of the optimization process.

Results

Double pipe

(A) (B) (C)

Double diffuser

The configuration of double diffuser case depicted in Figure 3.5, for which optimization process leads to different channels architecture when mixing the fluids and when separating them, is reconsidered using the modified inverse permeability coefficient of equation (3.8) to separate fluid subdomains. The same boundary conditions and meshing size of paragraph 3.2.2.2 are used in this example.

The result is similar to the optimization result using continuity objective function for fluid separation. However a slight difference still exists in the architecture of the channels, where power dissipation in the domain of Figure 3.12 is , whereas in Figure 3.7 the dissipation is equal to . In fact in the case of fluid separation using the continuity objective function, the method was able to preserve a more uniform width of both pipes and to get them closer by keeping only one solid mesh between fluid 1 and fluid 2. However both methods could be considered effective in optimizing channels shapes and separating fluid subdomains simultaneously. It should be noted also that both methods have the same performance regarding time calculation.

Constraint function

Implementation

The penalty objective function (3.1) is taken here as a constraint function instead of being inserted in the multi-objective optimization equation (paragraph 3.2). The optimization problem now is the following:

    1 1 2 1 2 1 2 1 2 1 2 2 , 3 min z( ), with z( ) (u, , ) subject to Equations (3.4) (3.5) (3.7) 1 ( , ) 1 1 ( , ) 1 G ( )
f f f f ij fc F f P T G d N Gd N f f d N                                   (3.15)
where and is the continuity function as defined in equation (3.1). and are respectively fluid 1 and fluid 2 maximum permissible volume; whereas is a very small number, i.e.

, to force continuity function to converge to 0, hence ensuring the separation of fluids subdomains. The color map is defined as follows: when the values go from 0 to 1, representative color goes from black, to dark red, then red, then yellow, then light yellow and finally to white. First Figure 3.14.A shows total convergence of field to a black and white solution; the two white pipes are made of 100 % fluid ( =1) and the remain cells are pure solid ( =0). Then it could be clearly seen from Figure 3.14.B, that full black color is present only in solid domain where there is no fluid at all, whereas for the pipe between the inlet and outlet section of fluid 1, the color is light yellow which means the fluid is not 100 % fluid 1 but there is a small portion of fluid 2 in the fluid mixture, otherwise the pipe color would had been white similarly to Figure 3.14.A. The remaining quantity of fluid 1 is clearly present in the upper pipe represented by dark red. Similarly in Figure 3.14.C the upper pipe that should correspond theoretically to fluid 2 only, is not made of pure fluid 2 and a small quantity of fluid 2 could be seen in the lower pipe. Thus the main conclusion of the above results is that fluid subdomains are not fully separated. In fact in fluid 1 field, the lower and upper pipes are made respectively of 95.5 % and 4.5 % of fluid 1, and vice versa for fluid 2 field. It could be concluded that the presence of continuity function as a constraint function in the optimization problem resulted in the convergence of field to an upper limit and a lower limit , whereas the same function taken as an objective function in paragraph 3.2 allowed to reach 0 and 1, thus total separation of fluids subdomains.

Results

The second drawback of the method is the extremely high number of iterations needed to reach the global optimal solution. In the example considered above, the algorithm reached the them near outlet region. Figure 3.16 shows topology optimization result in bi-fluid topology optimization method, with and without separation of fluids subdomains Continuity objective function method was used in the case of fluid phases separation (Figure 3.16.B). Without considering separation of fluid subdomains, the result of the optimization method is similar to the mono-fluid topology optimization result obtained in literature, as seen in Figure 3.16.A. In this case, all fluid cells are made of a single equivalent fluid, which is a mixture of 50 % fluid 1 and 50 % fluid 2 (

). However when continuity objective function is included in the multi-objective function, each fluid is transported in its own straight pipe, as seen in Figure 3.16.B. This structure has the lower pressure drop in case of forcing each fluid to enter and leave the domain from its predefined inlet and outlet boundaries.

Double pipe with AR=1.5: counter flow

The example of the previous paragraph is repeated, with two fluids flowing in opposite directions. The flowing direction of fluid 2 is inverted; it enters the domain from the east edge and leaves it from the west edge, as seen in Figure 3.17. All remaining geometrical characteristics, flow conditions and optimization parameters are similar to the previous paragraph. In Figure 3.19 fluid power dissipation throughout the optimization process is sketched for the cases of parallel and counter flows with fluid separation. The red curve shows a smooth and fast convergence to the optimal solution. However, the blue curve shows an extremely high number of oscillations between iterations 1 and 80. The conflict between objective functions leads to these oscillations: fluid power dissipation objective function forces the fluid to take the shortcut path between inlet and outlet boundaries at the same edge, whereas continuity objective functions builds an horizontal solid layer to forbid the fluid from taking the shortcut path. However after iteration 80, two horizontal pipes are created and the rest of iterations adjusts the solid-fluid interface in each pipe to finally converges to two parallel and straight pipes, as presented in Figure 3.18.B.

Conclusion

This chapter focused on the separation of fluid 1 and fluid 2 subdomains, an essential part that should be implemented in the algorithm presented in chapter 2, hence to be considered "bifluid" topology optimization method. Three different methods were implemented and tested. The first one by adding a continuity function to the multi-objective function that penalizes the existence of two different fluids in the same cell or neighbors cells. The second one by adding velocity absorption term to the inverse permeability coefficient of Darcy's law, which increases the power dissipation function when two different fluids are present in the same cell or in neighbor cells. Finally the last method consisted of taking the continuity function considered in the first method as a constraint function instead of an objective function. The first two methods showed approximately similar performance regarding the convergence of field, total separation of fluid subdomains, calculation time and forbidding the fusion of different fluid's pipes in the domain. However the first method was slightly better than the second one in the convergence of field (pipe architecture). The method of constraint function is found to be less attractive due three important drawbacks. First its inability to completely separate fluid subdomains, where a small fluid mixture will remain in fluid pipes. Second, the high number of iterations needed to reach the optimal solution. Finally, the method will begin to separate fluids after setting up the architecture of the pipes, which will lead in some configuration to join them instead of forcing each fluid to take its own independent path from the first stages of the optimization process. In the next chapter, maximization of heat transfer rate between the fluids will be considered, which will lead to the design of a fluid-fluid heat exchanger device optimized topologically using density method.

Chapitre 4: Transfert de masse et de chaleur en optimisation topologique bi-fluide

Introduction

Dans chaque échangeur de chaleur, un phénomène physique se produit toujours: pour chaque diminution infinitésimale de la différence de températures entre les deux fluides à la sortie de l'échangeur de chaleur, correspond une augmentation infinitésimale de la perte de charge quand la géométrie optimale est changée. Il existe donc toujours un compromis entre l'amélioration du transfert de chaleur et la chute de pression dans un échangeur de chaleur. De ce comportement se pose la question suivante: quelle est la structure qui assure une chute de pression minimale pour un taux de transfert de chaleur requis entre les fluides? Ou pour une puissance de pompage fixe, quelle est la structure qui fournit la différence minimale entre la température de sortie du flux chaud et celle du flux froid? Ces questions mettent en place un problème d'optimisation mathématique afin d'obtenir une réponse, qui est la géométrie optimale.

Fonction multi-objectif

La fonction objectif utilisée pour la maximisation du transfert de chaleur est celle qui correspond à la quantité récupérée par le fluide froid ou celle dissipée par le fluide chaud. Cependant pour des raisons numériques concernant l'analyse de sensibilité, le taux de transfert de chaleur sera calculé sur les deux fluides, donc la fonction objectif de récupération d'énergie exprimée en (4.3) sera égale au double de la quantité de chaleur transmise.

La solution d'optimisation multi-objectif qui correspond à la minimisation de la perte de charge et de maximisation de la récupération d'énergie, est nommée la solution optimale de Pareto. L'ensemble des solutions optimales de Pareto forment la frontière de Pareto, à partir de laquelle une solution devrait être choisie.

Il existe deux grandes familles d'optimisation multi objectif : l'articulation de préférence à priori qui consiste à choisir la préférence de chaque objectif avant le lancement du processus d'optimisation, et l'articulation de préférence à posteriori qui consiste à choisir une solution de la courbe de Pareto après l'optimisation. Dans ce travail, la méthode des sommes pondérées (4.4) qui appartiennent à la famille d'articulation de préférence a priori est adoptée. La partie critique de l'approche de la somme pondérée est le choix des poids où il pourrait être difficile de distinguer si l'effet des poids dans l'optimisation est de compenser les différences de grandeur des fonctions ou d'indiquer l'importance relative d'un objectif. Pour cela on a recours à la normalisation des fonctions objectifs, ainsi la fonction multi-objectif (4.6) unique à minimiser ne sera pas dominée par une fonction objectif au détriment d'une autre, ce qui permettra aux facteurs de pondération de refléter plus précisément les préférences. La dérivée totale de la fonction multi-objectif par rapport aux variables de design est exprimée dans (4.10). , l'optimiseur va changer la géométrie afin de diminuer la distance solide entre le fluide chaud et le fluide froid. Le tableau 4.6 montre l'augmentation de l'angle de courbure des canaux et la longueur de la section droite et la diminution de la distance entre les tuyaux quand diminue. Ces paramètres structurels témoignent de l'opposition entre la maximisation de perte de charge et la maximisation du transfert de chaleur pour les structures optimales. Pour , l'augmentation du poids de transfert de chaleur produit des géométries de plus en plus onduleuses, ce qui entraine des augmentations très importantes dans la fonction de dissipation de puissance dans le fluide. La Figure 4.24 montre la courbe de Pareto, qui correspond à l'ensemble des solutions optimales, c.à.d. la géométrie qui assure le maximum de transfert de chaleur pour une perte de chaleur fixe et vice-versa. In every heat exchanger device a physical phenomenon always occurs: for every infinitesimal decreasing in temperature difference between the two fluids at the exit of the heat exchanger, corresponds an infinitesimal increasing in pressure drop. Hence a trade off always exists between heat transfer enhancement and pressure drop in a heat exchanger. From this behavior arises the following question: what is structure that has the minimal possible pressure drop for a required heat transfer rate between the fluids? Or for a giving pumping power, what is the structure that provides the minimal temperature difference at outlet sections of hot and cold streams? These questions set up a mathematical optimization problem in order to be answered In this chapter, heat transfer will be considered with and without separation of fluid phases. First a single pipe will be studied without fluids separation to show the distribution of fluids having different thermal conductivities inside the pipe. Second, fluid phase separation will be considered in the maximization of heat transfer between the fluids, which will force each fluid to have his own channel.

Transfert de masse et de chaleur en optimisation topologique bi

Objective function 4.2.1. Heat transfer rate

For an incompressible fluid flowing through a solid domain having temperature at inlet and at exit, and not submitted to phase change, the total thermal power recovered is the following [START_REF] Yunus | Thermodynamics An Engineering Approach[END_REF]:

  21 P Q m h mC T T     (4.1)
Replacing by its expression: where is the area surface, the velocity and the density, equation (4.1) is written in its integral form as follows for :

    uu PP outlet inlet Q C T d C T d         (4.2)
The objective function representing maximization of heat transfer is the amount of heat recovered by the cold fluid or amount of heat dissipated by the hot fluid. However for numerical reasons regarding the sensitivity analysis, the heat transfer rate will be computed on both fluids, hence the objective function will be equal to the double of and is expressed as follows:

   

Multi-objective function

Fluid power dissipation and heat transfer objective functions present a trade-off, in which the optimization of the first one will be at the detriment of the second, and vice versa. This means any optimal solution design that leads to a decreasing in pressure drop, will lead to a decreasing in thermal energy transferred from hot fluid to cold fluid, and conversely. Hence the solution from optimizing simultaneously pressure drop and heat transfer in multi-objective optimization procedure is known as Pareto Optimal solution. Pareto Optimal is defined as follows in [118]: A vector is considered Pareto Optimal if there is not another solution vector that would minimize an objective function ( , without simultaneously causing an increase in at least one other objective function. Pareto is defined by the following statement: "Take from Peter to pay Paul" in order to show that the advantage in an objective function is taken from the pocket of the second objective function. The set of Pareto optimal solutions form the Pareto frontier, from which one solution should be picked. The choice of desirable point from Pareto curve is up to the decision maker. In literature there exist two main families of methods to choose a single solution [START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF]: Methods a priori articulation of preferences where a preference for every objective function is determined before running the optimization. Thus the final single solution reflects the preferences settled a priori. The second family of methods is a posteriori articulation of preferences that allows drawing the entire Pareto curve, than to choose manually a single solution from the curve. In this work a priori articulation of preference strategy is adopted. A widely used method in multi-objective optimization with a priori articulation of preferences is the weighted sum method defined as follows [START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF]:

11 1 nn i i i ii F w f w    (4.4)
where is the number of objective functions and are the weights that reflect the preference of every function.

The optimization problem aims now to minimize the single objective function known as aggregate objective function. The critical part of weighted sum approach is the choice of weights , where it could be difficult to distinguish if the effect of the weights in the optimization is to compensate the differences in objective-function magnitudes or to indicate the relative importance of an objective [START_REF] Marler | The weighted sum method for multi-objective optimization: new insights[END_REF]. Many methods were suggested in literature for assessing weight [START_REF] Tzeng | Multiple Attribute Decision Making: Methods and Applications[END_REF]. However the most used method is to transform the functions, in order to have the same magnitudes in all functions. Then the aggregate objective function can't be dominated by an objective function at the detriment of another which will allow the weights factors to reflect more accurately the preferences. The transformation procedure, known as normalization of objective functions, allows them to vary between zero and one. The normalization procedure is the following [START_REF] Grodzevich | Normalization and Other Topics in Multi-Objective Optimization[END_REF]:

min max min ii i ii ff f ff    (4.5)
where are the limits of Pareto optimal set and are determined as follows: minimization of alone gives and maximization of alone gives . Finally combining equations (4.4) and (4.5), the multi-objective function is expressed as follows in the case of fluid separation using continuity function: (1 )

d d e e cc d d e e f f f f F w w w f f f f f       (4.6)
It should be noted that does not need normalization because its bounded between 0 and 1. Furthermore, is not included in , because continuity objective function is a penalty function that has no contradiction with the other objectives when the problem converges to the optimal solution. Objective function of heat transfer has a minus sign because it is maximized, whereas fluid power dissipation and continuity functions are minimized.

Sensitivity analysis

The gradient of with respect to design variables is computed according to the method presented in paragraph 2.7, with the adjoint vector and residual function have the following forms:

U U U V V V P P P T T T AX B AX B J AX B AX B                                   (4.7)
The matrices and are expressed as follows:

V P JJ JJJ J J U V P J J J J J X U V P J J J JJ U V T                                                                            (4.8)
The matrices of partial derivatives of with respect to , and are computed according to paragraph 2.7.1.2.

In the calculation of , the velocity is located on the boundaries of the design domain but temperature needs to be interpolated between the values at fluid boundaries and the temperature values at the first cell inside the design domain in order to get the temperature exactly at the domain boundary. Hence have a partial derivative with respect only to the temperature field. The vectors and are expressed as follows:

0 0 0 0 0 ee e ff X f T                  (4.9)
Finally after computing the total derivative , the gradient of the multi-objective function is computed as follows: It should be noted that in the case of minimization of energy recovery, there is no trade off between the two objectives because and can be minimized simultaneously.

is a constant to control the preference of each objective, such that fluid distribution corresponds to the minimal fluid power dissipation, and affects only the distribution of the two fluids within the pipe. Finally, the equivalent thermal conductivity in each design cell and all other physical parameters are computed using equation (2.19). The fluid power dissipation, which is similar for the three cases, is . Energy recovery for case A is , for case B is and for case C is . It can be observed from the result of case B, that low conductive fluid is not distributed uniformly along the interface of the pipe and the central core. In fact, fluid 1 phase is thicker near west and east edges of the central core, than at near its north and south edges. This explains that the horizontal heat flux inside the central core has a much higher negative effect on fluid energy recovery than the vertical heat flux. The vertical heat flux in the central core takes In this case, the algorithm used the low conductive fluid to split the central core into two parts, to the same reason to limit the effect of the horizontal heat flux discussed above.

126

A study was conducted in order to show the conflict between using a fluid quantity to create an insulation layer in the middle of the central core (case C), or use the same quantity in the main flow regions to decrease pressure drop and drive fluid pipes toward heat source (case B). The latter also increases , but less than the increment caused by splitting the central core. The study consist of computing energy recovery objective function for (100 values of were considered) in optimal structures of cases B and C, than to compute the total multi objective function for

The assumption made in this study is to consider a fixed upper bound of energy recovery for all values of . The results are sketched in Figure 4.12. It can be seen from the curves that for , optimal structure of case B is not sufficient anymore to provide the required heat transfer rate, and optimal result of case C has a lower objective function despite having a much higher fluid power dissipation.

u i f i f T , 2 , 2 u i f i f T , 1 , 1 u o f o f T , 2 , 2 u o f o f T 127 and
. Maximum porosities of both fluids are . Solid thermal conductivity is . The problem will be solved first by minimizing fluid power dissipation only, and then maximization of heat transfer will be considered with and without separation of fluids subdomains, which will lead to different architectures.

Minimization of fluid power dissipation

Heat and mass transfer without separation of fluids subdomains

In this study, equation (4.6) is used as the multi-objective function to be minimized, with (without fluid separation). Topology optimization results are presented in Figure 4.15. The results are characterized by two main features. First the creation of two parallel pipes between fluid 1 inlet and fluid 2 outlet sections and between fluid 1 inlet and fluid 2 outlet sections. Second, when decreases, the pipes are driven toward south and north edges in order to increase solid thickness between the pipes, which means minimization of heat transfer rate. To explain this behavior, we consider, , , , and are respectively average temperatures at fluid 1 inlet, fluid 1 outlet, fluid 2 inlet and fluid 2 outlet sections, and average inlet and outlet velocities at both fluids inlet and outlet boundaries. Heat transfer objective function (equation (4.3)) could be approximated as follows:

    , 1 , 1 , 2 , 2 e p o f i f p o f i f f C U T T C U T T      (4.13)
For and , equation (4.13) becomes:

    , 1 , 2 , 1 , 2 273 423 150
e p o f o f e p o f o f f C U T T f C U T T          (4.14)
Hence, increases when increases and decreases. Normally should correspond to the outlet temperature of the cold fluid ( ) and should correspond to the outlet temperature of the hot fluid (

). However the distribution of inlet and outlet boundaries of both fluids, allows the creation of continuous pipes between inlet and outlet sections of different fluids, without intersection of the pipes. Thus, there is a great advantage for the maximization of , that hot fluid leaves the domain from outlet boundary of fluid 1, and the cold fluid leaves the domain from the outlet boundary of fluid 2. Furthermore, in this case increasing and decreasing requires minimization of heat transfer rate, thus in case of the two pipes are driven away from each other, as seen in Figure 4. [START_REF] Rao | Effect of flow distribution to the channels on the thermal performance of a plate heat exchanger[END_REF] However, it will be seen later that in fluid-fluid heat exchangers a small increasing in heat transfer objective function requires a huge increasing in fluid power dissipation. It should be noted that in this example has no physical meaning and doesn't represent heat transfer rate between the streams. In fact, optimization results in Figure 4.15 are a clear example of failure of the optimization problem of heat and mass transfer in bi-fluid topology optimization in absence of continuity objective function, which ensures that each fluid enters and leaves the domain from its predefined boundaries without any intersection with the other fluid. Hence fluid-fluid heat exchangers cannot be designed and optimized using mono-fluid topology optimization method, even for the same fluids properties. . For iterations 1 to 105, the two vertical pipes on east and west edges are driven closer to each other to increase heat transfer rate. Until this point represents heat transfer load from the hot stream to the cold stream. Between iterations 105 and 110, the two vertical pipes merge together, and both objective functions undergoes a fast increasing in their values. After iteration 110, the fluids are transported by two horizontal pipes, which are driven away from each other between iterations 110 and 314. This increases by a small quantity because it does not represent heat transfer rate anymore and is only valid from a mathematical point of view. It also leads to a huge decreasing in pressure drop because the pipes become straight.

Heat and mass transfer with separation of fluids subdomains

In this paragraph fluid subdomains separation will be included in the optimization problem of heat and mass transfer, using continuity objective function method ( in equation (4.6)). The results of topology optimization for are sketched in Figure 4.17. We can notice that a single pipe for each fluid is generated, connecting its inlet and outlet boundaries located on the same edge of the domain. Using continuity objective function ensured the convergence to a realistic physical solution, in which each fluid takes its own path independently from the second, hence represents the heat transfer rate between the streams. The variations of and throughout the optimization process are sketched in Figure 4.18. The first initial configuration is the optimal solution of minimization of fluid power dissipation alone (Figure 4.14.B). The graph shows that increases progressively until reaching a trade-off solution between and for . It can be noticed that for every infinitesimal increasing in , corresponds an infinitesimal increasing in , which is a normal and realistic physical behavior in heat and mass transfer, as discussed in the first chapter. It should be noted that and have increased respectively by 1791% and 175% in relation to the optimal configuration in case of minimization of pressure drop alone (Figure 4.14.B). This shows that pressure drop and heat transfer don't vary proportionally; in fact a small decreasing in outlet temperature difference between the fluids requires a much higher increasing in pumping power. decreases, the thickness of solid wall separating the fluids is reduced by driving the pipes closer to each other. At the closet distance between the pipes, a single solid cell separates fluid 1 and fluid 2 cells. For the residence time of fluid in the domain is increased, by increasing the total distance of the pipe in which the fluid is transported. The optimal structure of , is similar to the structure of with lower diameter and higher total length, in order to increase furthermore heat transfer rate. To maximize heat transfer between fluids, both pipes should have the minimal solid thickness for the longest distance possible . Decreasing and increasing requires a higher bending pipe angle . Higher bending angle leads to higher pressure drop, hence higher fluid power dissipation. Thus, the trade-off between geometrical parameters reflects the trade off between heat transfer maximization and pressure drop minimization. To show the effect of bending pipe on pressure drop, fluid power dissipation function is computed between sections A and B for pipes of optimal structures of and , as seen in However for this structure, the increasing in heat transfer is very small in relation to the much higher increment in pressure drop.

A) (B) (C) (D) (E) (F) (G) (H) (I) (J) (K) (L)
A) (B) (C) (D)
structures is bending angle at inlet sections, between the fluid pipe and a horizontal line normal to inlet boundary. This angle is equal to 26° for the optimal structure of , 35° for the optimal structure of and 47° for the optimal structure of .

 For , heat transfer is further maximized by increasing heat transfer surface and generating structures that have similar effect as fins in heat exchangers. When decreases, cold and hot pipes became longer and thinner which minimize the outlet temperatures difference but also causes a huge increasing in pressure drop. It could be seen from Figure 4.26, that wavy shapes are first generated near inlet boundaries of hot and cold fluids, where the fluids have their maximum and minimum temperatures respectively, than those wavy shapes are intensified throughout the whole domain when decreases. In the chart's legend fd and fe stand respectively for fluid power dissipation and energy recovery, and pf and cf correspond respectively for parallel flow and counter flow. The graph shows that fluid power dissipation curves for both fluid arrangements have approximately similar values. However heat transfer is higher for counter flow especially for low values of . This shows clearly that counter flow arrangement is able to provide higher heat transfer rate than parallel flow arrangement, with similar and sometime The graph shows that for the entire range of between 0 and 1, the outlet temperatures of cold stream are higher in counter flow and the outlet temperatures of hot streams are higher in parallel flow, which means a higher heat transfer rate in counter flow. The difference between outlet temperatures in parallel and counter flows ( and ) increases when decreases. . In parallel flow, fluid pipe has two bending sections, near outlet and near inlet regions, to bring both pipes closer together in the center of the domain. Both pipes in the center of the domain are horizontal. In counter flow, fluid pipes have bending sections only in near inlet region, and after the bending section the fluid leaves the domain by a straight but tilted pipe.

To test the performance of both structures in parallel and counter flows, we considered the optimal structure of in parallel flow sketched in Figure 4.31.A and denoted as structure A, and the optimal structure of in counter flow sketched in Figure 4.31.B and denoted as structure B. These structures were considered for comparison because they have similar values of , for structure A and for structure B. In parallel flow, for and , in structure A and for structure B. Thus was only higher by 0.06%, which could be rather caused by numerical error, hence both structures could be considered with same affectivity in parallel flow. In counter flow, for same inlet temperatures cited above, in structure A, and in structure B. However despite the difference is only by 0.6% between both structures, but is still more considerable than in case of parallel flow, and was detectible by the optimization algorithm, which explains the convergence to structure B in counter flow.

Conclusion

In this chapter, maximization of heat transfer at the detriment of fluid power dissipation in bi-fluid topology optimization was considered. This has allowed design and optimization of fluid-fluid heat exchangers using advanced mathematical optimization technique of topology optimization. First, distribution of fluid 1 and fluid 2 within the same pipe was studied, in minimization and maximization of heat recuperation by the entire fluid channel crossing the domain, whose is subjected to a constant temperature at its south and north edges. The results showed the validation of the proposed bi-fluid method in simultaneous optimization of heat and mass transfer without fluids separation. It should be noted that these results are valid from a theoretical point of view, and their application in real life engineering devices may not be possible even with immiscible fluids. Then, heat transfer maximization between separated fluid subdomains was studied. Fluid 1 and fluid 2 subdomains formed respectively the cold and hot streams. Continuity objective function ensured that each fluid is transported in its own pipe that connects its predefined inlet and outlet boundaries. An example of possible failure structure in absence of fluid phases separation was presented. Results have shown that separation of fluids subdomains is essential in bi-fluid topology optimization, not only for the determination of the type of fluid in each pipe, but also to ensure a convergence to a realistic engineering solution. Finally, optimization of double pipe heat exchangers in parallel and counter flows was studied. The minimization of the multi-objective function allowed drawing the optimal Pareto curve, which represents the trade-off between minimization of pressure drop and maximization of heat transfer in design and optimization of heat exchangers. Enfin, une nouvelle fonction objectif qui calcule le taux de chaleur transféré entre les fluides était ajoutée au problème d'optimisation multi-objectif pour optimiser le transfert de chaleur et la perte de charge simultanément. Premièrement, le problème était résolu sans séparation des fluides pour étudier le positionnement des deux fluides ayant des conductivités thermiques différentes pour la minimisation et la maximisation de la récupération d'énergie respectivement. Enfin pour les domaines comprenant des entrées et sorties différentes des deux fluides, la séparation des fluides était respectée dans le problème d'optimisation multi-objectif grâce à l'utilisation de la fonction de continuité. Les structures générées par l'optimiseur dépendent du paramètre de pondération qui détermine la préférence entre la maximisation de transfert de chaleur et la minimisation de la perte de charge.

Cependant, malgré le succès de la méthode proposée pour traiter le problème posé au début de la thèse, et pour établir la base d'une méthode de conception des échangeurs de chaleur spécifiquement, et divers dispositifs d'ingénierie qui inclut généralement le transfert de masse et chaleur, la méthode présente certaines limitations considérables exprimées comme suit :

 Le nombre extrêmement élevé des variables d'optimisation, qui est le double du nombre de variables de conception pour la même configuration dans les problèmes mono-fluide, ce qui limite l'application de la méthode sur des domaines relativement petits, ou forcent à utiliser des grands maillages. Ainsi la précision de la structure optimale finale diminue ainsi que la précision de la solution de la méthode des volumes finis.  La convergence de la méthode d'optimisation est fortement dépendante de plusieurs paramètres mathématiques dont on cite les plus importants : le paramètre de pénalisation de la fonction d'interpolation, les rayons des filtres de sensibilités et de conception, le poids de la fonction de continuité et les paramètres de la méthode des asymptotes mobiles.

Enfin les prochains travaux doivent étendre la méthode au domaine tridimensionnel et aux écoulements turbulents afin que la méthode soit utilisée par les concepteurs des échangeurs de chaleur dans des applications d'ingénierie réelles. Les contraintes de fabrications additives, et des méthodes numériques telles que l'adjoint continu et la GCMMA peuvent être également intégrées dans l'optimiseur. Chapter 5

Conclusion and Perspectives

Conclusion and limitations

The present thesis aims to extend the application of topology optimization on fluid flow problems to include two fluid subdomains in addition to a solid subdomain, whereas literature review showed that this kind of problems solved so far included only a single fluid in addition to the solid phase. To deal with the presence of two fluid phases and a solid phase in the optimization domain two interpolation functions were developed: mono-eta interpolation function that uses one design variable in each design cell based on normal distribution function, and bi-eta interpolation function that uses two design variables in each design cell, one for separating solid phase from both fluids phases and a second one for separating fluid 1 from fluid 2 phase in the fluid portion determined by the first design variable. Results showed that monoeta formulation depends strongly on initial guess whereas bi-eta formulation is capable of rearranging solid and both fluids freely without any limitations.

Various methods were then implemented to ensure the separation of fluids subdomains: penalty equation as an objective function, the same equation as a constraint function and using a physically unrealistic term added to the inverse permeability coefficient in the flow equation, which increases pressure drop when two fluids mix in the same channel. The penalty equation as an objective function and the modified inverse permeability coefficient methods showed similar efficiency in separating fluid subdomains and preventing them from mixing in the same cell or same channel, also they prevented the fusion or crossing of fluid channels where a minimal solid thickness is preserved between them. However the method of using penalty function as a constraint function was found to be not so efficient where it failed to completely separate fluids subdomains and was very time consuming due to the high number of iterations needed.

which will decrease the accuracy of CFD calculations and optimal structures results. In addition, calculation of Navier-Stokes and energy equations in 3D porous flows requires additional computational time in comparison with 2D CFD time calculation. In summary, moving from 2D to 3D domains in bi-fluid topology optimization introduces a large field of investigations and challenges in the development of the method due to the increasing of complexity and time calculation of direct CFD solver, discrete adjoint method and numerical optimization method.

Turbulent flow

All flow considered in this thesis are limited to laminar flows only. Above certain Reynolds number, flows become unstable where a chaotic and random state of motion is developed, called turbulent flow, observed in most of engineering practices. Contrary to laminar flow which is smooth and well organized and ordered, turbulent flow is unsteady and rotational flow structures are observed, where velocity and other flow properties vary randomly. As a result of those rotational structures, heat, mass and momentum are very effectively transfered in comparison with laminar flow. Thus despite the increasing in energy losses in turbulent flows, sometimes they are preferred over laminar flow in heat exchangers for their higher convective heat transfer coefficient. However application of topology optimization on turbulent flow requires the implementation and derivation of special equations and turbulence models, which are more complicated than CFD calculation of laminar flow due to the chaotic behavior of turbulent flows. Considering turbulent flows in topology optimization of 2D and 3D bi-fluid domains increases the complexity of the problem much more than moving from 2D to 3D in laminar flow (paragraph 5.2.1), and opens up a large field of investigation and studies, which could lead to much more promising and interesting works.

Numerical methods

Density method used in this thesis as topology optimization technique for heat and mass transfer in bi-fluid domain, is based on three main parts: a direct solver to solve Navier-Stokes and energy equations, a suitable method to compute the gradient of objective function with respect to design variables, and finally numerical optimization method that computes the optimal value of design variables by minimizing an objective function on the basis of gradient information. Finite Volume method (FVM), discrete adjoint method, and method of moving asymptotes (MMA) are respectively used in the three main parts of the algorithm. However the following methods should also be considered in the future:

 There are two methods to compute the gradient of the objective function using the adjoint technique: discrete adjoint method and continuous adjoint method. As already explained in paragraph 2.7, the main difference between the two methods is that in discrete method the physical equations are discretized and then derived whereas in continuous method they are derived then the derivative equations are discretized to solve the adjoint problem. The methods give close results especially with fine structured meshes, hence there is no preference of one of the methods in literature despite discrete adjoint method have better agreement with finite difference method. In topology optimization both methods are used in most of the works done in this field of study. In this thesis only discrete adjoint method is considered, thus continuous adjoint method should also be implemented. Then a comparison should be made between results, convergence process and calculation time for both methods.  Similarly to most of works done in the field of topology optimization, the Method of Moving Asymptotes (MMA) is used for numerical optimization. Globally Convergent Method of Moving Asymptotes (GCMMA), which is an upgraded version of MMA that ensures the convergence to global minimum, looks also appealing and has been used in topology optimization despite higher number of iterations and calculation time it needs. Hence the implementation of GCMMA in current design methodology developed in this thesis could be considered in the future.

Boundary conditions

The works conducted in this thesis were limited to fluid flow entering and leaving the optimization domain through a normal vector to the surface. However in real life engineering application, nothing can guarantee that fluid leaves the domain orthogonally which affect the shape of the optimal structure, hence the convergence procedure and the final optimal structure.

Manufacturing constraints

In first chapter it was stated that the advancement in additive manufacturing technology, pushed researchers and engineers to investigate and develop design methodologies using topology optimization techniques. In fact topology optimization may leads to complex structures unable to be manufactured using conventional techniques of material extraction and welding, whereas 3D printers solved nearly most of all problems related to manufacturing considerations and the limitations imposed on the field of topology optimization. However additive manufacturing still presents some constraints that have to be implemented in the optimization algorithm to be taken into consideration in the design process, otherwise the development of any efficient engineering device is useless if its manufacturability conditions are not respected. The manufacturing constraints in 3D printers are based on law of gravity applied on printing process, in which a piece cannot be printed in midair. In fact every layer of material printed serves as a support for the next layer during printing cycle. The following features that should be considered when designing and optimizing are the following:

 The minimum supported wall thickness: a constraint must be considered on the minimal walls, fins and other pieces inside the exchanger. This constraint depends on the minimal layer thickness the printer is able to print at each cycle (normally this thickness is around ).  Maximum overhang length for parts that includes an unsupported overhang, there is a maximum length for the overhang that could be printed without using a support, which range usually between and .  Maximum unsupported overhang angle without using a support, which is a variant constraint that depends on the material used. However parts with angles up to could be successfully printed.

Annex A1: Logarithmic mean temperature difference method

Heat exchangers sizing and rating are solved using Logarithmic mean temperature difference method based on the two following equation for heat flux transferred between the hot and cold fluids: where the heat flux, the temperature, the heat capacity which is the mass flow rate multiplied by the specific heat capacity, the overall heat transfer coefficient and the effective heat transfer surface is the logarithmic mean temperature difference that depends on flow arrangement and is expressed as follows:

  For counter flow arrangement:

For parallel flow arrangement:

1 , , 2 , , 1 , , 2 , ,

h i c o h o c i h i c i h o c o T T T T T T T T T T T T

                      (A1.4)
For all types of flow arrangements could be compute as follows:

,0 LM LM T F T    (A1.5)
where being the logarithmic mean temperature computed for counter flow arrangement and a correction factor corresponding to a specific flow type. is found in abacus of heat exchanger design books, or computed using explicit functions that depends in two parameters and defined as follows:

Annex A2: Effectiveness-Number of transfer unit method

One way of evaluating the performance of a heat exchanger is to compare it with an ideal, unrealistic heat exchanger where the temperature of the cold fluid at its outlet section reaches the inlet temperature of hot fluid. This could be only realized with a counter flow heat exchanger of infinite length. Hence the corresponding maximal heat transferred between the fluids is: For we will have and :

    max , ,
        , , , , max min , , min , ,

c c o c i h h i h o h i c i h i

c i C T T C T T q q C T T C T T       (A2.3)
is a dimensionless parameter that describes the performance of the heat exchanger and is expressed as following: Depending on the flow arrangement, a mathematical expression relating to and the ratio of lowest to highest heat capacity, can be developed (for detailed mathematical development of the functions see [START_REF] Shah | Fundamental of heat exchanger deisng[END_REF]). This expression is used for sizing and rating problems in heat exchanger design and does not need an iterative procedure. Examples for some common flow types: Parallel flow heat exchanger: Expressions for a wide range of heat exchangers technologies are summarized in [START_REF] Shah | Fundamental of heat exchanger deisng[END_REF].

Annex A3: Temperature effectiveness -number of transfer unit method

The heat transfer rate from cold to hot fluid is expressed as follows: was first developed for design of shell and tube heat exchangers before the introduction of , and was then adapted to different types of heat exchangers technologies.

Annex B1: Entropy generation method

Thermodynamic irreversibility is identified by computing the entropy generation due to heat transfer in the heat exchanger. The entropy generation is influenced by the temperature and the temperature difference distributions in the exchanger. The local temperature difference is the driving force for heat transfer in a heat exchanger and hence it influences the exchanger effectiveness. The entropy generation for an adiabatic open system is the sum of entropy generation of both fluids: where is the mass flow rate and is the entropy rate change between the inlet and outlet. Subscript 1 and 2 stand respectively for each of the hot stream and cold stream. By integrating between the inlet and outlet entropy generation due to heat transfer in heat exchanger becomes [START_REF] Shah | Fundamental of heat exchanger deisng[END_REF]: where the subscripts and stand respectively for inlet and outlet.

For an incompressible fluid the entropy generation caused by fluid friction is [START_REF] Shah | Fundamental of heat exchanger deisng[END_REF]: where

    , , , , , , , , , ln ln
The problem often encountered with the shape optimization based on the minimization of entropy generation number is the "entropy generation paradox" where for some flow conditions and boundary conditions the decreasing in entropy generation number yield to a reduction in heat exchanger effectiveness in contrary to expected result. Entropy generation paradox was clearly demonstrated by Shah et Skiepsko [START_REF] Shah | Entropy Generation Extrema and Their Relationship With Heat Exchanger Effectiveness-Number of Transfer Unit Behavior for Complex Flow Arrangements[END_REF] by analyzing the relationship between heat exchanger efficiency and entropy generation in 18 heat exchangers with different structures, the results showed that when the entropy generation reaches the extremum, the efficiency of the heat exchangers can be at its maximum, minimum or anything between.

To avoid entropy generation paradox Guo et al. [START_REF] Guo | Optimization design of shell-and-tube heat exchanger by entropy generation minimization and genetic algorithm[END_REF] have introduced the modified entropy generation number by using the ratio of heat flow and inlet temperature of cold fluid to nondimensionalization of entropy generation. The entropy generation number becomes [START_REF] Guo | Optimization design of shell-and-tube heat exchanger by entropy generation minimization and genetic algorithm[END_REF]:

.

gen cold inlet s ST N Q   (B1.6)
The results showed that for the optimization problems where the heat load is given both the initial entropy generation number (B1.5) and the newly introduced one (B1.6) lead to the same results. However if the duty heat is not fixed entropy generation number defined by Bejan (B1.5) suffered from the entropy generation paradox, which is avoided by the new formulation (B1.6).

Annex B2: Entransy dissipation method

Based on analogy between heat flux and electrical current, Guo et al. [START_REF] Guo | Entransy-A physical quantity describing heat transfer ability[END_REF] have developed in 2007 a new physical quantity called entransy to describe heat transport potential capacity of an object. The name entransy came from "en-transy" where "en" stands for energy and "transy" stands for transport. The definition of entransy G is:

1 2 G UT  (B2.1)
By analogy to electricity, the entransy of an object is the heat transfer ability which corresponds to the electrical energy in a capacitor that describes its charge transfer ability. the internal energy in equation (1.8) corresponds to electrical quantity in a capacitor and the temperature corresponds to the voltage. The entransy is transported during heat transfer similarly to the transport of electric energy during electric conduction. By this definition it is obvious to show that the heat transfer ability of an object depends on its potential to transfer heat which is the temperature and its capacity to be able to transfer it which is .

Guo et al. have developed the expression of entransy dissipation for heat conduction rate defined as [START_REF] Guo | Entransy-A physical quantity describing heat transfer ability[END_REF]:

. h qT     (B2.2)
Where is the heat flow density and the temperature. In a heat transfer process the thermal energy is conserved but the entransy is not conserved like the entropy is generated, thus the entransy dissipation rate is a physical quantity for the measuring the irreversibility in heat transfer, thus used as objective function in shape optimization of heat exchangers.

For the heat conduction the entransy dissipation is computed by applying Fourier law to equation (B2.2) and integrating between inlet and outlet. The final form of entransy dissipation caused by heat conduction can be expressed as follows [START_REF] Guo | Entransy dissipation number and its application to heat exchanger performance evaluation[END_REF]: where the subscripts h,c,i and o stands respectively for hot, cold, inlet and outlet.

    2 2 2 2 , , , , 11 ( ) ( ) 22
The entransy dissipation related to fluid friction for an incompressible fluid is developed by applying the thermodynamic entransy dissipation expression for heat convection where a second term related to fluid viscosity appears in equation (B2.2), thus the entransy dissipation related to fluid friction in heat exchanger can be expressed as follows [START_REF] Xu | Application of entransy dissipation theory in heat convection[END_REF]: After replacing the cell faces properties values corresponding to QUICK scheme of equation D1.2 into the discretized equation of finite volume method, it appears that the main coefficient are not guaranteed to be positive, which is a condition to the interpolation scheme to be stable and ensures the convergence of the iterative solver. It's found that under certain conditions, stability problems may occur with QUICK scheme and unbounded solutions may appear. The maximum value of Peclet number that ensures stability is 8/3, thus QUICK scheme as presented by Leonard [126] is conditionally stable. However many authors have reformulated Leonard's scheme to alleviate stability problems by placing the terms that give rise to negative coefficients in the source term to retain all coefficients of transport properties of all nodes in the discretized equation positive. One of the best known QUICK reformulation is the one developed by Hayase et al [START_REF] Hayase | A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures[END_REF].The new formulation is achieved on the basis of satisfying five rules that guarantee physically realistic solution of the equations of conservation of mass and momentum approximated by the finite volume method. The discretized form of physical equations in Hayase et al QUICK scheme is the following:

    , , , , , , , , ln ln
                                                                
P P W W E E S S N N A A A A A S           (D1.3)
where the coefficients of equation D1. Using this approach, coefficients of equation D1.3 are always positive which satisfy the requirements of differencing schemes properties. At th iteration, the source term is evaluated using the values of physical properties of th iteration, hence it's differed by one iteration. However, after a large number of iterations both versions of QUICK schemes will converge to the same solution.

A A A A A C C C C A D C A D C A D C A D C S C C CC CC                                                                       

Mots Clés

Optimisation topologique, méthode de densité, bi-fluide, écoulement laminaire, échangeurs de chaleur Abstract Heat exchangers are devices widely used in various energy systems. The present research work focuses on the design of single-phase bi-fluid heat exchangers by using topology optimization methods. Unlike conventional size and shape optimization methods, topology optimization methods allow greater design freedom and do not require prior definition of the exchanger geometry.

Hence, bi-fluid topology optimization consists of freely reorganizing two fluids and one solid in the optimization domain. Both fluids should connect inlet sections to outlet sections while avoiding any fluid mixture inside the domain. SIMP method "Solid Isotropic Material with Penalization" is used within the framework of this thesis. This method includes various algorithms that were formulated and tested: finite volume method was selected for solving the direct physical problem, discrete adjoint method was used for the calculation of the gradient of the objective function, and the method of moving asymptotes was adopted to guide the numerical optimization. Simulation results have allowed the definition of various heat exchanger shapes in 2D. The number of inlet and outlet as well as the flow rates of each fluid have been varied. The works have shown the ability of this method to design innovative shapes of heat exchangers. Hence, the thesis establishes the basis of a new design methodology of heat exchangers.

Figure 1 . 1 :

 11 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 1 :

 21 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 2

 22 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 3

 23 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 .

 2 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 6 :

 26 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 7 :

 27 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 8 :

 28 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 9

 29 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 12 :

 212 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 13 :

 213 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 14 :

 214 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 15 :

 215 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 16 :

 216 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 17 :

 217 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 .

 2 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 .

 2 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 .

 2 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 22 :

 222 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 . 23 :

 223 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

Figure 2 .Figure 2 . 25 : 88 Figure 3 . 2 : 90 Figure 3 . 3 : 92 Figure 3 . 5 : 92 Figure 3 . 6 : 93 Figure 3 . 7 : 93 Figure 3 . 8 : 94 Figure 3 . 9 : 96 Figure 3 . 10 : 97 Figure 3 . 11 : 102 Figure 3 . 16 : 103 Figure 3 . 17 :Figure 4 . 14 :

 222588329033923592369337933894399631097311102316103317414 Figure 1.1: Fluid to fluid heat exchanger .................................................................................................................. Figure 1.2:Examples of heat exchangers construction technologies [6] ..................................................... Figure 1.3: Common vortex generators used [16] ............................................................................................ Figure 1.4: General algorithm for heat exchanger optimization ................................................................. Figure 1.5: Comparison of size, shape and topology optimization ............................................................. Figure 1.6:Size, shape and topology optimization applied on heat exchangers ................................... Figure 1.7: Pareto optimal front in shape optimization[39] ......................................................................... Figure 1.8: Example of resulting blades from tube tank heat exchanger shape optimization[39]19 Figure 1.9: Topology optimization on mechanical structures problems [47] ....................................... Figure 1.10: Discretized domain in topology optimization. .......................................................................... Figure 1.11: Topology optimization in heat and mass transfer ................................................................... Figure 1.12: Topology optimization results in Stokes flow [66] .................................................................. Figure 1.13: Design of a bend pipe using level set method [78] ................................................................... Figure 1.14: Design of a diffuser using level set method [78] .......................................................................Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem .....................Figure 2.2 Main Loop .....................................................................................................................................................Figure 2.3 Inner Loop ....................................................................................................................................................Figure 2.4 Normal Distribution function for and . ........................................................ Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient ...............................Figure 2.6: Constraint functions for fluid 1 and fluid 2. ..................................................................................Figure 2.7: Interpolation function (2.10) for different values of .............................................................Figure 2.8:Bi-eta interpolation function ................................................................................................................Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell ......................................... Figure 2.10: Representation of a control volume ....................................................................................... Figure 2.11: Distribution of transport quanitites and their coefficients .................................................. Figure 2.12: Initial configuration of double pipe example for mono-eta formulation .......................Figure 2.13: Topology optimization results for various initial values of .............................................Figure 2.14: Mono-eta interpolation function for double pipe configuration .......................................Figure 2.15: Bi-eta topology optimization of double pipe domain .............................................................Figure 2.16: Variation of power dissipation objective function throughout the optimization process .................................................................................................................................................................................Figure 2.17: Variation of power dissipation objective function from iterations 24 to 60 .........Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and (B) ...........................................................................Figure 2.19: Topology optimization result for ............................................................................Figure 2.20 Topology optimization of single pipe with .......................................................... Figure 2.21: Initial configuration of diffuser case .............................................................................................Figure 2.22: Calculation of global power dissipation .......................................................................................Figure 2.23: Calculation of local power dissipation ..........................................................................................Figure 2.24 Topology optimization results of the diffuser case ..................................................................

  4: Thermal and hydraulic performance of topology optimization results for various values of weighting factor .................................................................................................................................... 131 Table 4.5: Thermal and hydraulic performance of topology optimization results for various values of ...................................................................................................................................................................... 134 Table 4.6: Geometrical parameters and of for structures of .............................. 137 Table 4.7: Thermal and hydraulic performance of topology optimization results for various values of weighting factor .................................................................................................................................... 139 xi Nomenclature Heat flow Mass flow

L

  'optimisation topologique dans les problèmes d'écoulement était initialement limitée à de faibles nombres de Reynolds et à un état stationnaire (écoulement de Stokes) sans tenir compte des effets d'inertie. Ensuite, divers auteurs ont étendu la procédure d'optimisation pour couvrir une plus large gamme de nombre de Reynolds, des effets d'inertie (flux de Darcy-Stokes et de Navier-Stokes), des forces corporelles non uniformes et des flux instationnaires.

Figure 1 . 1 :

 11 Figure 1.1: Fluid to fluid heat exchanger

Figure 1 . 2 :

 12 Figure 1.2:Examples of heat exchangers construction technologies [6]

Figure 1 . 3 :

 13 Figure 1.3: Common vortex generators used[START_REF] Jacobi | Heat transfer surface enhancement through the use of longitudinal vortices: A review of recent progress[END_REF]

Figure 1 .

 1 Figure 1.4: General algorithm for heat exchanger optimization

Figure 1 . 5 :

 15 Figure 1.5: Comparison of size, shape and topology optimization

Figure 1 . 6 :

 16 Figure 1.6: Size, shape and topology optimization applied on heat exchangers

Figure 1 . 7 .

 17 Examples of resulting blades solutions from the Pareto front are sketched in Figure 1.8.

Figure 1 . 9 :Figure 1 . 10 :

 19110 Figure1.9: Topology optimization on mechanical structures problems[START_REF]Topology Optimization Program[END_REF]

Figure 1 .

 1 Figure 1.11 represents typical topology optimization results in heat and mass transfer problems. Figure 1.11.A represents an application of topology optimization on pure heat conduction problems. The final solution shows the optimal distribution of limited amount of high conductive material (in black) to evacuate the heat generated in low conductive material (in white) through a small patch (heat sink) located in the middle of the left boundary of the domain[START_REF] Marck | Topology optimization using the SIMP method for multiobjective conductive problems[END_REF].

Figure 1 . 11 :

 111 Figure 1.11: Topology optimization in heat and mass transfer

3 .

 3 Coupling heat transfer to mass transfer problem in case of two fluids, which correspond to the design of a fluid to fluid heat exchanger.

Figure 1 . 12 :

 112 Figure 1.12: Topology optimization results in Stokes flow[START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF]

Figure 1 . 14 :

 114 Figure1.13: Design of a bend pipe using level set method[START_REF] Duan | Shape-topology optimization for Navier-Stokes problem using variational level set method[END_REF]

  b) Fonction d'interpolation bi-etaDans [66] Borrvall et Peterson ont introduit la première fonction d'interpolation aux problèmes d'écoulement des fluides. Ils ont également fourni les preuves mathématiques que la fonction permet la convergence vers une solution optimale. De nombreux auteurs ont adopté ce schéma d'interpolation dans lequel correspond à la phase solide, et à la phase fluide et à un milieu poreux ayant des paramètres physiques intermédiaires entre le fluide et le solide, comme le montre l'équation (2.15). La convexité de la courbe de la fonction est ajustée à l'aide d'un paramètre de pénalisation qui joue un rôle décisif dans la convergence de la méthode d'optimisation (Figure 2.7). D'autre part, Sigmund [52] a proposé un schéma d'interpolation pour les problèmes de structure mécaniques multi-matériaux (2 matériaux et le vide). Ce schéma consiste à interpoler d'abord entre les deux matériaux, puis entre le matériau équivalent et le vide. En appliquant la même logique, le schéma d'interpolation bi eta consistera alors à interpoler premièrement entre le fluide 1 et le fluide 2 et ensuite interpoler entre le fluide équivalent et le solide, comme le montre l'équation (2.18). L'expression introduite par Borrvall et Peterson pour les problèmes mono fluide est remplacé dans l'expression générique du schéma bi-fluide (2.18). Ainsi deux variables de conception et définissent la nature de la phase dans une maille : pour le solide, pour le fluide 1 et pour le fluide 2. La forme finale de la fonction d'interpolation bi-eta bi-fluide est présentée dans l'équation(2.19). Enfin les fonctions contraintes de porosité qui limitent les porosités maximales du fluide 1 et fluide 2 sont formulées en fonction de et , comme le montre l'équation(2.21)

Figure 2 .

 2 15 montre la convergence de à son optimum global. Cependant dans toutes les mailles fluides est égale à 0.5, ce qui signifie que chaque maille fluide consiste en un mélange 50% fluide 1 et 50% fluide 2. En fait, ce résultat était attendu puisque les deux fluides ont des propriétés physiques et des limites de porosités égales, ce qui signifie que la séparation des fluides n'au aucun effet sur la fonction objectif ou contraintes. Dans la deuxième étude, un déséquilibre est créé dans les propriétés et les limites de porosités des deux fluides afin d'observer la convergence de . La configuration initiale est représentée dans la Figure 2.18.A. Trois cas de calcul ont été considérés :  Cas A : , et . des cas A, B et C sont respectivement présentées dans Figure 2.18.B, Figure 2.19 et Figure 2.20. Pour le cas A, a convergé vers 0 et 1 alors que e est égale à 0.5 dans tout le domaine vers 0 et 1, ce qui signifie que la convergence des mailles fluides vers le fluide 1 ou le fluide 2. Le minimum de perte de charge correspond à l'accumulation du fluide 1 au milieu du canal d'écoulement et le fluide 2 dans le reste du canal. Dans le cas C, pour pouvoir utiliser toutes les quantités du fluide 1 et du fluide 2, prend des valeurs entre 0.71 et 0.74 dans toutes les mailles fluides. Ainsi les valeurs de était proches du rapport de limite de porosité du fluide 2 et la porosité totale des deux fluides (fluide 1 + fluide 2). Ces trois exemples montrent la possibilité de la convergence des champs et simultanément, donc la possibilité de distribuer de façon optimale deux fluides et un solide dans le domaine d'optimisation en utilisant la fonction d'interpolation pénalisée bi-eta.

  Walls at constant temperature .  Walls subject to normal heat flux .  Adiabatic walls .  Inlet flows with parabolic profile at constant temperature for first fluid.  Outlet flows with parabolic profile u 01 .  Inlet flows with parabolic profile at constant temperature for second fluid.  Outlet flows with parabolic profile .

  Figure 2.1: Initial guess and final solution of bi-fluid topology optimization problem

Figure 2 . 2

 22 Figure 2.2 Main Loop

Figure 2 . 3

 23 Figure 2.3 Inner Loop

Figure 2 .

 2 Figure 2.4 Normal Distribution function for and .

Figure 2 . 5

 25 Figure 2.5 Mono-eta interpolation function for inverse permeability coefficient

Figure 2 . 6 :

 26 Figure 2.6: Constraint functions for fluid 1 and fluid 2.

Figure 2 . 8 :

 28 Figure 2.8:Bi-eta interpolation function

Figure 2 . 9

 29 Figure 2.9 Constraint function for fluid 1 (A) and fluid 2 (B) in a single cell

Figure 2 . 11 : 57 Figure 2

 211572 Figure 2.11: Distribution of transport quanitites and their coefficients

Figure 2 . 12 :

 212 Figure 2.12: Initial configuration of double pipe example for mono-eta formulation

Figure 2 .

 2 [START_REF] Jang | The thermal-hydraulic characteristics of staggered circular finned-tube heat exchangers under dry and dehumidifying conditions[END_REF].B and C showed that topology optimization results reached in cases B and C are not global optimums, and their corresponding values of the objective function is higher than in case D represented in Figure2.13.D, which is the global optimal solution for this problem.

Figure

  Figure

Figure 2 Table 2 . 1 :

 221 Figure 2.13.A 0.00081224 Figure 2.13.B 0.00030124 Figure 2.13.C 0.00030125 Figure 2.13.D 0.00021481

Figure 2 . 13 :Figure 2 . 14 :

 213214 Figure 2.13: Topology optimization results for various initial values of

  optimization algorithm was able to easily reach the global optimal solution and to distribute fluid 1 and fluid 2 simultaneously but with the help of the initial guess.

1 .Figure 2 . 15 :

 1215 Figure 2.15: Bi-eta topology optimization of double pipe domain

Figure 2 . 16 :

 216 Figure 2.16: Variation of power dissipation objective function throughout the optimization process

Figure 2 .Figure 2 .Figure 2 .

 222 Figure 2.18 Single pipe configuration (A) and its optimisation result (B) ( field) using bi-eta interpolation function for and

Figure 2 .Figure 2 .

 22 Figure 2.19 show the optimization results in case where , and . The presence of only white and black color, and the absence of any grey cells in Figure 2.19.B and Figure 2.19.C show the complete separation of fluid 1 and fluid 2 phases, which means complete converge of field to 0-1 solution. Fluid 1 with lowest permeability among the three phases is accumulated in the middle of the fluid channel and fluid 2 which has is accumulated in the remaining space between fluid 1 and solid-fluid interface at the upper and lower boarders of the channel.Figure 2.19.D shows the three phases simultaneously in the

Figure 2 .

 2 [START_REF] Caputo | Heat Exchanger Optimized Design Compared With Installed Industrial Solutions[END_REF].A shows the complete convergence of field to black-and-white solution, where the shape of the fluid channel is similar to the channel in the structure of Figure 2.18.B. Figure 2.20.B and Figure 2.20.C show that the fluid inside the channel is a mixture of fluid 1 and fluid 2 with the major quantity belongs to fluid 2. In fact to take advantage of the excess porosity of fluid 2, varies between 0.71 and 0.74 in all fluid cells which have . This means that in average each cell is made of 72.5 % of fluid 2 and 27.5 % of fluid 1, which is approximately equal to the ratio of fluid 2 mean porosity over the calculation domain:

  the solid phase. The Reynolds number based on inlet section length and the average velocity at the inlet section is . The domain has 96 mm side length and is discretized into design cells.

Figure 2 . 22 : 78 Figure 2 . 23 :

 22278223 Figure 2.22: Calculation of global power dissipation

Figure 2 . 25 :

 225 Figure 2.24 Topology optimization results of the diffuser case

Figure 2 .

 2 Figure 2.26 Topology optimization results of bend pipe case

Figure 2 Table 2 . 2 :

 222 Figure 2.26.A Figure 2.26.B Figure 2.26.C

Chapitre 3 :

 3 Séparations des sous-domaines fluides 1. Introduction La méthode d'optimisation topologique bi-fluide présentée dans le 2 ème chapitre manque encore d'une caractéristique très importante: sa capacité à séparer les deux fluides. Par conséquent la méthode ne peut pas encore être considérée comme une méthode d'optimisation de topologie bi-fluide capable de créer des canaux indépendants de chaque fluide. Les résultats du chapitre précédent ont montré que le problème bi-fluide est encore résolu comme dans le cas d'optimisation topologique mono-fluide, avec un seul fluide équivalent, dans le cas où les propriétés physiques et les limites de porosité du fluide 1 et fluide 2 sont égales. L'objectif de ce chapitre est de trouver la méthode la plus appropriée pour :  Premièrement, empêcher le mélange du fluide 1 et du fluide 2 dans la même cellule lorsque l'algorithme converge vers la solution finale.  Deuxièmement, assurer la continuité de chaque fluide, c.à.d. que chaque fluide soit transporté dans un ou plusieurs tuyaux qui relient les sections d'entrées et de sorties de ce fluide, qui sont prédéfinies par le concepteur sur les limites du domaine.  Troisièmement, s'assurer que les canaux de deux fluides différents ne puissent pas se mélanger ou se croiser, et que l'algorithme doit être capable de maintenir une épaisseur solide minimale entre les canaux.

  En partant d'un domaine initial homogène d'un matériau intermédiaire, la méthode d'optimisation topologique bi-fluide doit séparer les deux fluides et le solide en minimisant la fonction objectif et en respectant les fonctions contraintes. Le calcul de coefficient inverse de perméabilité par la fonction d'interpolation pénalisée assure la convergence de à 0 ou 1, c.à.d. à la phase solide ou fluide, sans distinction entre fluide 1 et fluide 2. Cependant il n'y a aucune raison mathématique ou physique à faire converger vers 0 ou 1, c.à.d faire converger le fluide vers le fluide 1 ou fluide 2, ce qui fait que toutes les cellules fluides sont constituées d'un mélange de 50% fluide 1 et de 50% fluide 2. Par conséquent, une fonction objectif est formulée et ajoutée au problème d'optimisation multi-objectif afin d'empêcher mathématiquement le mélange des fluides. La fonction objectif, appelée fonction de continuité est exprimé dans les équations (3.1),(3.2),(3.3)et (3.4). La Table

  .2.A représente le domaine d'optimisation du premier cas d'étude ayant une entrée et une sortie pour chaque fluide. La solution finale, représentée dans la Figure 3.2.B, montre la création de deux canaux parallèles : le premier reliant l'entrée et la sortie du fluide 1 formé par des mailles 100% fluide 1 ( et ) et un deuxième reliant l'entrée et la sortie du fluide 2 formé par des mailles 100% fluide 2 ( et ). Ainsi l'architecture des canaux générés est similaire à celle du paragraphe 2.10.2.1, et la seule différence demeure dans la nature du fluide dans chaque canal. Ainsi, la fonction de continuité avait un effet uniquement sur les variables de conception , c.à.d sur la composition des mailles fluides, et n'avait aucun effet sur , c.à.d sur la forme des canaux d'écoulement. Cependant, ce n'est pas le cas pour toutes les configurations, où pour certaines géométries une contradiction entre la fonction de dissipation de puissance du fluide et celle de la continuité peut se produire. Cela conduira à différentes architectures de canaux de fluide avec et sans séparation des fluides, comme le montre le deuxième cas d'étude. La deuxième configuration est présentée dans la Figure 3.5. Les résultats d'optimisation topologique avec et sans séparation des fluides par la fonction de continuité, sont présentés respectivement dans les Figure 3.6 et 3.7. Sans séparation des fluides, l'algorithme joint les canaux d'écoulement en un seul canal au milieu du domaine, puis le sépare avant la sortie. Or, en considérant la fonction de continuité dans la fonction multi-objectif, chaque fluide est transporté par son propre canal reliant l'entrée et la sortie de ce fluide, indépendamment du canal d'écoulement de l'autre fluide. Au moins, une maille solide sépare les canaux d'écoulement des différents fluides. Cet exemple montre que l'utilisation de la fonction objectif de continuité assure la séparation totale des fluides et omettra toutes les valeurs intermédiaires de . L'algorithme recherchera la forme des canaux de fluide pour avoir la perte de charge minimale mais en respectant la condition de séparation des fluides.

  Le domaine d'optimisation étudié est rectangulaire ayant un rapport d'aspect AR=1.5. Les porosités maximales des fluides sont , ce qui correspond à un tuyau droit rempli de fluide entre l'entrée et la sortie de chaque fluide. Le problème est résolu avec et sans séparation des fluides pour des écoulements à co-courant et à contre courant. Les configurations initiales des deux arrangements d'écoulement sont représentées dans les Figure 3.15 et 3.17. Le résultat d'optimisation topologique du cas à co-courant montre que sans séparation des fluides, l'optimiseur a joint les deux canaux dans la zone voisine des sections d'entrées des deux fluides et les a séparé dans la zone voisine des sections de leurs sorties (Figure 3.16.A). Cependant pour le même cas d'écoulement mais avec séparation, l'optimiseur assure le transport de chaque fluide par son propre canal en créant deux canaux parallèles, comme le montre la Figure 3.16.B. Pour l'écoulement à contre courant sans séparation des fluides, deux raccourcis reliant les entrées et les sorties des différents fluides sont crées, respectivement dans les régions voisines des frontières gauche et droite (Figure 3.18.A). Enfin, la Figure 3.18.B montre que la fonction de continuité évite la création de ces raccourcis, ainsi la séparation des fluides est assurée dans l'écoulement contre courant. La minimisation de la dissipation de puissance avec séparation des fluides en écoulement parallèle a nécessité 59 itérations, alors qu'à contre-courant il a fallu 374 itérations pour mettre en place la solution finale. Le nombre élevé d'itérations à contre-courant reflète la difficulté de l'algorithme d'optimisation à éviter le raccourci du canal du fluide, car les entrées et les sorties sont du même côté et proches l'une de l'autre. La Figure 3.19 montre les fonctions de dissipation de puissance tout au long du processus d'optimisation pour les cas des flux à co-courants et à contre courants avec séparation des fluides. La courbe rouge, celle du cas à co-courant, montre une convergence facile et rapide vers la solution optimale finale. Cependant la courbe bleue, celle du cas à contre courant, montre un nombre extrêmement élevé d'oscillations entre les itérations 1 et 80. Le conflit entre les fonctions objectifs conduit à ces oscillations: la fonction objectif de dissipation de puissance du fluide oblige le fluide à prendre le raccourci entre les bornes d'entrée et de sortie, alors que la fonction objectif de continuité construit une couche solide horizontale pour interdire au fluide de prendre le chemin du raccourci. Cependant, après l'itération 80, deux tuyaux horizontaux sont créés et le reste des itérations ajuste l'interface solide-fluide dans chaque tuyau pour finalement converger vers deux tuyaux parallèles et droits.

Figure 3 . 2 :

 32 Figure 3.2: Double pipe example with fluid separation using continuity objective function

Figure 3 . 3 : 92 Figure 3 Figure 3 . 5 :

 3392335 Figure 3.3: Variation of power dissipation objective function throughout the optimization process

  presented in Figure3.5 has a rectangular shape of width and length and a maximum porosity of 0.2 for each fluid. The domain is discretized into 100 x 150 square meshes of 1mm side length. The inlet areas for each fluid are equal and three times bigger than their respective outlet sections, hence the velocity at the exit is u u u u to ensure the mass flow continuity through the domain. The Reynolds number at the inlet is .

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Results of double diffuser example without continuity objective function

Figure 3 . 8 :

 38 Figure 3.8: Summary of separation of fluids subdomains using the inverse permeability coefficient

14 )Figure 3 .Figure 3 . 9 :

 14339 Figure 3.8 summarizes the method of fluids subdomains separation using the inverse permeability coefficient.

Figure 3 . 10 :

 310 Figure 3.10: Comparison of fluid dissipation function for different fluid mixtures

Figure 3 .

 3 Figure 3.10 compares the values of power dissipation function between the global optimal solution (Figure 3.10.A) and other undesired solutions. In Figure 3.10.B the fluids are permuted, where fluid 1 flows through the inlet and outlet sections of fluid 2 initially predefined, and vice versa for fluid 2. The fluid 1-fluid 2 interfaces at the inlets and outlets boundaries of fluids have increased by 335 %. For the same reason, is increased by 692 % in Figure 3.10.C when inserting a small portion of fluid 1 in the channel of fluid 2 and a small portion of fluid 2 in the channel of fluid 1, which will act like an obstacle to the flow. Finally in Figure 3.10.D, the fluid flowing in both pipes is made of 50 % fluid 1 and 50 % fluid 2 (). It could be clearly seen that the fluid power dissipation function increases dramatically by 37517 %. This explains why the algorithm will begin to separate the fluids from the first stages of the optimization process by converging the values of to 1 near the inlet and outlet sections of fluid 2 and to 0 near the inlets and outlets sections of fluid 1, and later along the entire pipes as long as they are created progressively by the convergence of field. However small convergence steps in MMA method and a density filter on field are applied to avoid the convergence to local optimal solution similar to Figure3.10.C.

Figure 3 . 11 :

 311 Figure 3.11: Comparison of power dissipation due to fluid friction, with and without fluid separation

Figure 3 .

 3 Figure 3.11 compares the minimization of in cases of using the initial permeability coefficient () used in double pipe example in paragraph 2.10 and the actual inverse permeability coefficient used in this example to separate fluids subdomains. It should be noted that in both cases, the domain is initialized with the same fields of and . For the initial configuration, in case of using the inverse permeability coefficient of equation (3.8), whereas when the physical problem is solved with the initial

Figure 3 . 12 :

 312 Figure 3.12: Topology optimization result for double diffuser configuration with the modified inverse permeability coefficient

Figure 3 . 13 :Figure 3 .

 3133 Figure 3.13: Initial domain configuration

Figure 3 . 16 :

 316 Figure 3.16: Topology optimization results of double pipe example in parallel flow, (A) without fluids phases separation and (B) with fluid phases separation.

Figure 3 .

 3 Figure 3.18 shows the topology optimization results with and without fluid separation. In case of no fluid separation, it is clearly seen that two pipes are created, respectively near west and east domain edges. The first pipe, begins at fluid 1 inlet boundary and ends at fluid 2 outlet boundary, both on the west edge. Similarly, the second pipe begins at fluid 2 inlet boundary and ends at fluid 1 outlet boundary, both on east edge. Both pipes are made of the equivalent fluid (50 % fluid 1 and 50 % fluid 2). Joining inlet and outlet sections of different fluids was permitted

  -fluide sans séparation des fluides a) Minimisation de récupération d'énergie Dans ce paragraphe la séparation des fluides n'est pas prise en compte afin d'étudier la capacité de la méthode à répondre au bon positionnement des fluides ayant des conductivités thermiques différentes dans un seul canal d'écoulement. La configuration comprenant une seule entrée et une seule sortie pour les fluides est représentée dans la Figure 4.1. La température du fluide à l'entrée est Les frontières nord et sud du domaine sont soumis à une température fixe Les propriétés du fluide 1, du fluide 2 et du solide sont résumées dans le tableau 4.1. Pour minimiser la récupération de chaleur par le fluide durant son écoulement à travers le domaine, l'algorithme préférera le fluide ayant la conductivité thermique la plus faible, donc une quantité plus limitée de fluide 2 est autorisée à suivre sa distribution optimale. Ainsi les limites maximales de porosité sont et . Dans le cas de minimisation de récupération d'énergie, le signe de est inversé, ainsi la fonction multiobjectif dans ce cas est exprimée dans (4.11) Les résultats de l'optimisation topologique sont présentés dans la Figure 4.2. Le fluide 2 ayant une conductivité thermique 100 fois inférieure à celle du fluide 1, est accumulé sur les périphéries inférieures et supérieures du tuyau. En fait, le fluide 2 agit comme une couche isolante pour empêcher la chaleur d'atteindre les régions à haute vitesse au milieu du canal. La valeur de de ma solution finale est comparée avec celle des deux autres cas : premièrement en inversant la conductivité des deux fluides en laissant le même positionnement et deuxièmement en prenant , c.à.d comme si on avait un seul fluide ayant une conductivité thermique équivalente . Dans les deux cas, la récupération d'énergie a augmenté, ce qui montre que l'optimiseur était capable de trouver le meilleur positionnement des deux fluides dans le tuyau. b) Maximisation de récupération d'énergie La même configuration est reconsidérée afin de maximiser la récupération d'énergie, ainsi la fonction multi objectif minimisée est celle exprimée en (4.6) avec (sans séparation). Les trois cas d'études sont les suivants : d'entrée du fluide est et celle des parois nord et sud est Dans cette étude le problème d'optimisation comprend une seule contrainte, celle de la porosité minimale du solide exprimée en (4.12), avec . Les Figure 4.5 et 4.6 montrent les résultats de l'optimisation topologique pour les cas A et B respectivement. Pour le cas A, l'égalité des propriétés du fluide, dans toutes les cellules fluides signifie que chaque cellule est constituée d'un fluide équivalent constitué d'un mélange de 50% de fluide 1 et 50% de fluide 2 ( ). Cependant, dans le cas B l'optimiseur a utilisé les deux fluides ayant des conductivités thermiques différentes ( et ) sans avoir aucun mélange des deux fluides dans la même cellule de conception. Pour les deux cas le tuyau d'écoulement est divisé en deux petits tuyaux pour se rapprocher de la source de température et augmenter la récupération d'énergie (champs similaire). Cependant pour le cas B, la distribution des deux fluides n'était pas homogène comme pour le cas A. Le fluide 1 a la faible conductivité thermique pour isoler thermiquement le noyau central afin que la chaleur récupérée par le fluide ne soit pas dissipée dans le noyau central pour augmenter sa température. La Figure 4.9 montre la diminution de la fonction de récupération d'énergie quand la conductivité thermique de fluide isolant augmente. La Figure 4.10 montre le résultat d'optimisation topologique du cas C. La solution optimale finale montre qu'une barrière de fluide est créée au milieu du noyau central afin d'empêcher le refroidissement du fluide chaud à la sortie par le fluide froid à l'entrée à travers ce noyau ayant une haute conductivité thermique. La Figure 4.11.A montre qu'il n y a pas un écoulement dans la barrière centrale, ainsi cette quantité de fluide n'a aucun effet sur la diminution de fonction de dissipation de la puissance du fluide. Par contre la Figure 4.11.B montre l'importance de cette barrière centrale afin d'éviter le flux de chaleur horizontale inverse à travers le noyau central. Les solutions finales des cas B et C montrent une contradiction entre l'utilisation d'une quantité du fluide pour créer une barrière centrale ou utiliser cette quantité dans les régions principales de l'écoulement pour diminuer la perte de charge et rapprocher le fluide de la source de température. La variation de la fonction multi objectif pour en fonction de la conductivité thermique du fluide isolant (Figure 4.12) montre que la solution du cas B est meilleure seulement pour . 4. Transfert de masse et de chaleur en optimisation topologique bi-fluide avec séparation des fluides La configuration initiale pour ce paragraphe est présentée dans la Figure 4.13. Les fluides 1 et 2 sont considérés respectivement comme étant les débits froid et chaud. Les températures d'entrées des deux fluides sont respectivement et . Le reste des propriétés physiques des deux fluides est identique. Les porosités maximales des deux fluides sont . Le problème est résolu premièrement sans séparation des fluides ( ). Les résultats de l'optimisation dans la Figure 4.15 montrent la création de deux tuyaux parallèles : le premier reliant l'entrée du fluide 1 à la sortie du fluide 2 et le deuxième reliant l'entrée du fluide 2 à la sortie du fluide 1. En effet, afin d'augmenter , la température de la sortie flux chaud doit diminuer et celle du fluide froid doit augmenter, c.à.d. maximiser la température de sortie du fluide 1et minimiser la température de la sortie du fluide 2 . Cependant la distribution des sections d'entrée et de sortie des deux fluides permet la création de deux tuyaux parallèles entre les frontières des différents fluides sans aucune intersection entre ces tuyaux. Ainsi, afin de maximiser , le fluide chaud quitte le domaine de la section de sortie du fluide 1, et le fluide froid quitte le domaine de la section de sortie du fluide 2. Dans ce cas, l'augmentation de et la diminution de nécessitent la minimisation du taux de transfert de chaleur. Ainsi pour les deux tuyaux sont éloignés l'un de l'autre.Dans ce qui suit, la fonction de continuité est utilisée pour assurer la séparation des fluides ( ). La Figure4.17 montre la solution finale pour . Il est clair qu'un tuyau est généré dans la région gauche pour le fluide froid (fluide 1) et dans la région droite pour le fluide chaud (fluide 2). L'utilisation de la fonction objectif de continuité a assuré la convergence vers une solution physique réelle, dans laquelle chaque fluide prend son propre chemin indépendamment de l'autre, ainsi représente la quantité de chaleur transférée du débit chaud au débit froid.Les résultats d'optimisation topologique pour plusieurs valeurs de sont présentés dans la Figure 4.19. Les fonctions objectifs et et les températures de sortie du fluide 1 et du fluide 2 des structures optimales de la Figure 4.19, sont présentées dans le tableau 4.4. Lorsque w diminue, l'épaisseur de la paroi solide séparant les fluides est réduite en rapprochant les tuyaux les uns des autres. Pour , la longueur des canaux est augmentée et leur diamètre est diminué afin de maximiser le transfert de chaleur. 5. Cas d'étude : double tuyaux Dans cet exemple, le fluide 1 et le fluide 2 entrent et sortent du domaine par leurs sections d'entrée et de sortie prédéfinies à la limite du domaine. La séparation des fluides dans le domaine est assurée en utilisant la fonction objectif de continuité. Par conséquent, l'équation (4.6) est utilisée comme fonction multi-objectif unique à minimiser. Deux arrangements d'écoulement ont été considérés: écoulement à co-courant dans lequel les fluides chauds et froids entrent au niveau de la frontière gauche et quittent le domaine au niveau de la frontière droite, et l'écoulement à contre courant en inversant le sens d'écoulement du fluide chaud. Toutes les propriétés des deux fluides sont égales, ainsi que leurs limites de porosité . La température d'entrée du fluide froid (fluide 1) est et celle du fluide chaud (fluide 2) est . a) Ecoulement à co-courant : La configuration initiale est présentée dans la Figure 4.20. Les résultats d'optimisation topologique pour plusieurs valeurs de poids de pondération sont présentés dans la Figure 4.21 et le tableau 4.5. Pour

  b) Ecoulement à contre-courant :La configuration initiale est présentée dans la Figure4.[START_REF]Equipartition of entropy production. An optimality criterion for transfer and separation processes[END_REF]. Les résultats d'optimisation topologique pour plusieurs valeurs de poids de pondération sont présentés dans la Figure4.26 et le tableau 4.7. Les résultats sont semblables au cas de l'écoulement à co-courant : en donnant un poids plus important au transfert de chaleur au détriment de la perte de charge, l'optimiseur va rapprocher les tuyaux des différents fluides au début, puis la création des ondulations pour les petites valeurs de , en diminuant le diamètre des canaux et en augmentant la longueur des tuyaux dans le domaine. La Figure 4.27 montre la courbe de Pareto pour le cas à contre courant. c) Comparaison entre l'écoulement à co-courant et à contre courant : La Figure 4.28 montre la variation des fonctions objectifs de perte de charge et de récupération d'énergie en fonction de pour les géométries optimales dans les deux cas d'arrangement des écoulements. Le graphique montre que les courbes de dissipation de puissance de fluide pour les deux arrangements ont des valeurs approximativement similaires. Cependant, le transfert de chaleur est plus important pour le cas à contre-courant, en particulier pour les faibles valeurs de . Ceci montre clairement qu'un agencement à contre-courant est capable de fournir un flux de chaleur entre les fluides plus élevé que l'agencement à écoulement parallèle, avec une chute de pression similaire et parfois inférieure. Par conséquent, un écoulement à contre-courant est préféré à un écoulement parallèle dans un échangeur de chaleur fluide-fluide. La variation de la température moyenne aux sections de sortie du fluide froid (fluide 1) et du fluide chaud (fluide 2) en fonction du facteur de pondération , pour les écoulements parallèles et à contre-courant, est illustrée dans la Figure 4.29. Le graphique montre que pour toute la plage de entre 0 et 1, les températures de sortie du flux froid sont plus élevées à contrecourant et les températures de sortie des flux chauds sont supérieures en flux parallèle, ce qui signifie un taux de transfert de chaleur plus important pour l'écoulement à contre courant.

  vector normal to the surface, which means is always positive at fluid outlet boundaries and negative at fluid inlet boundaries. and denotes for fluid inlet and outlet sections, respectively for fluid 1 and fluid 2. It should be noted that fluid 1 is always considered the cold fluid, and fluid 2 the hot fluid to have a positive value of .

1 .

 1 Single pipe: Minimization of energy recovery

Figure 4 . 1 : 1 :

 411 Figure 4.1: Configurations of single pipe case

  Because both objectives are minimized, the multi objective function in this case is the following:

Figure 4 . 2 :

 42 Figure 4.2: Topology optimization result for single pipe by minimization of

Figure 4 . 3 :

 43 Figure 4.3: Temperature and velocity profiles at velocity exit boundary for cases A, B and C for R e =10

Figure 4 . 4 :

 44 Figure 4.4: Temperature and velocity profiles at velocity exit boundary for cases A, B and C for R e =500

Figure 4 . 5 : 2 : 122 Figure 4 Figure 4 . 6 :

 452122446 Figure 4.5: Topology optimization result of single pipe by maximization of : case A

  Figure 4.9 shows the decreasing of thermal energy recovery function with the increasing of the thermal conductivity of fluid 1.

123 placeFigure 4 . 7 :Figure 4 . 8 :Figure 4 . 9 :Figure 4 . 11 :

 123474849411 Figure 4.7: Temperature fields of optimal results of cases A and B

  Figure 4.10 shows the topology optimization result of case C, which is similar to case B with

  Figure 4.11.A shows that there is no velocity in the vertical fluid barrier in the central core. This means that this quantity of fluid used has no effect on decreasing fluid power dissipation, and has only effect on the temperature field, as seen in Figure 4.11.B. The fluid power dissipation function in this case is , which is higher by 48%, in relation to fluid power dissipation of the optimal result of case B. However heat recovery is , which means an increment by 25.9 % in relation to case B.

Figure 4 . 12 :

 412 Figure 4.12: Variation of in function of for optimal structures of cases B and C for

4. 5 .

 5 Heat and mass transfer in bi-fluid topology optimization with fluids separation 4.5.1. Initial configuration

Figure 4 . 13 :

 413 Figure 4.13: Initial configuration of double pipe with fluid inlet and outlet boundaries on the same edge

Figure 4 . 14 :

 414 Figure 4.14: Topology optimization results: Minimization of fluid power dissipation without fluid separation (A) and with fluid separation (B)

Figure 4 . 15 :

 415 Figure 4.15: Topology optimization results of heat and mass transfer without fluid separation for (A) and (B)

Figure 4 . 16 :

 416 Figure 4.16: Variation of and throughout the optimization process for without fluid separation

Figure 4 . 17 :

 417 Figure 4.17: Topology optimization result of heat and mass transfer with fluid separation for

Figure 4

 4 Figure 4.19: Topology optimization result of heat and mass transfer with fluid separation for

Figure 4 . 21 :Figure 4 . 22 :Figure 4 . 23 :

 421422423 Figure 4.21: Double pipe configuration with parallel flow arrangement: structure variation with respect to the weighting factor

  Figure 4.23. For the structure of sketched in Figure 4.23.A, fluid power dissipation in bending section is l t 

Figure 4 . 24 :

 424 Figure 4.24: Double pipe with parallel flow arrangement: Pareto frontier

Figure 4 . 25 :

 425 Figure 4.25: Initial configuration of double pipe with counter flow arrangement

  Figure 4.26 shows combined subdomains of solid, fluid 1 and fluid 2 for the optimal structures that correspond to different values of weighting factor . Power dissipation and heat transfer objective functions, hot and cold fluids outlet temperatures and temperature difference for the structures depicted in Figure 4.21 are summarized in Table 4.7.

Figure 4 . 27 : 6 . 3 .

 42763 Figure 4.27: Double pipe with parallel flow arrangement: Pareto frontier

Figure 4 .

 4 Figure 4.28 summarizes fluid power dissipation and heat transfer functions of pareto frontiers for parallel and counter flow. In the chart's legend fd and fe stand respectively for fluid power dissipation and energy recovery, and pf and cf correspond respectively for parallel flow and counter flow. The graph shows that fluid power dissipation curves for both fluid arrangements have approximately similar values. However heat transfer is higher for counter flow especially for low values of . This shows clearly that counter flow arrangement is able to provide higher heat transfer rate than parallel flow arrangement, with similar and sometime

  x10 -7 W) lower pressure drop. Hence counter flow is preferred over parallel flow in fluid-fluid heat exchanger.

Figure 4 . 28 :

 428 Figure 4.28: Fluid power dissipation and heat transfer objective function variation in function of weighting factor for parallel and counter flows arrangements

Figure 4 . 29 :

 429 Figure 4.29: Outlet temperatures of cold and hot streams variation in function of weighting factor for parallel and counter flows arrangements

Figure 4 .

 4 Figure 4.29.The graph shows that for the entire range of between 0 and 1, the outlet temperatures of cold stream are higher in counter flow and the outlet temperatures of hot streams are higher in parallel flow, which means a higher heat transfer rate in counter flow. The difference between outlet temperatures in parallel and counter flows ( and ) increases when decreases.

Figure 4 . 30 :

 430 Figure 4.30: Outlet temperature of cold and hot streams throughout optimization process in parallel and counter flows for

Figure 4 .Figure 4 .

 44 Figure 4.30 shows the variation of average temperature at outlet boundaries of hot and cold fluids throughout the optimization process of in parallel and counter flows. The fluid power dissipation in optimal structures of both cases are respectively and . The graph shows that throughout the majority of iterations, counter flow provides a lower temperature difference between average outlet temperatures of cold and hot streams. The temperature difference between the streams are 73.98°C and 54.63°C respectively at the final iteration of parallel and counter flow cases.

Chapitre 5 :

 5 Conclusion et perspectivesDans cette thèse, l'optimisation topologique dans les problèmes de transfert de masse et de chaleur a été étendue pour des domaines bi-fluides, optimisant ainsi la distribution de 2 fluides et un solide dans le domaine. Deux fonctions d'interpolation pour forcer la matière fictive intermédiaire dans chaque maille à converger vers l'une des 3 phases (fluide1, fluide 2 ou solide) ont été formulées et testées. La première fonction mono-eta utilise une seule variable de conception dans chaque maille et est basée sur la superposition de trois fonctions de distribution normales. La seconde, bi-eta, utilise deux variables de conception dans chaque maille, pour interpoler entre le fluide 1 et le fluide 2 afin de calculer la propriété physique du fluide équivalent, et pour interpoler entre le fluide équivalent et le solide. Ainsi définit la forme des canaux d'écoulement et définit la nature du fluide dans chaque maille de conception à l'intérieur des tuyaux. La deuxième formulation a été trouvée plus efficace et a été adoptée dans le reste de la thèse en raison de la dépendance de la première vis-à-vis de l'estimation initiale. Les résultats ont montré la capacité de la formulation bi-eta à répondre au bon positionnement des fluides pour la minimisation de la fonction de dissipation de puissance dans le fluide.Ensuite l'optimiseur avait besoin d'assurer le transport de chaque fluide dans le domaine par son propre canal indépendamment de l'autre fluide. Ce dernier devrait connecter les entrées et les sorties du fluide correspondant. La séparation des fluides devrait également assurer qu'une couche solide d'épaisseur minimale d'une maille de conception sépare les canaux d'écoulement des différents fluides. Trois méthodes étaient implémentées dans l'algorithme d'optimisation afin de pénaliser la présence des fluides différents dans la même maille ou dans des mailles voisines. Dans la première méthode, la séparation était assurée par l'intermédiaire d'une fonction de continuité prise comme fonction objectif. Dans la deuxième le coefficient inverse de perméabilité est modifié afin de séparer les fluides en considérant un effet mur entre les différents fluides. Et enfin dans la troisième méthode la fonction de continuité était prise comme fonction contrainte. Les deux premières méthodes ont montré une efficacité similaire en séparant les deux fluides et en les empêchant de se mélanger dans la même cellule ou le même canal. Cependant la troisième méthode a été jugée moins efficace puisqu'elle ne permettait pas de séparer complètement les fluides et nécessitait beaucoup de temps en raison du nombre élevé d'itérations nécessaires.

  the temperature of the fluid flowing outside the tubes and for the fluids inside the tubes.

5 )

 5 effectiveness and is the difference between inlet temperatures of hot and cold fluids. Temperature effectiveness is a dimensionless parameter computed for each fluid in function of its corresponding number, heat capacity rate ratio and the flow arrangement:Finally sizing and rating problem in design of heat exchangers are solved by computing unknown variables by using the relation for the corresponding flow arrangement and definitions of dimensionless parameters in A3.3.

3 )

 3 Finally the total entropy generation in the heat exchanger is:Actually the objective function used to optimize heat exchangers is the entropy generation number introduced by Bejan [123] defined as follows:

Figure A :

 A Figure A: Representation of the control volume and its neighbor nodes

Table 2 .

 2 1: Objective function of topology optimization result of Figure 2.13 ....................................... 68 Table 2.2: Investigation of topology optimization results of the bend pipe case ................................. 80 Table 3.1: Continuity function coefficients at extremum values of and ....................................... 89 Table 4.1: Properties of fluid 1, fluid 2 and solid ............................................................................................ 117 Table 4.2: Investigation of topology optimization results .......................................................................... 121 Table 4.3: Thermal and hydraulic performance of topology optimization results for and .

............................................................................................................................................................................. 128

Table 4 .

 4

Table 1 .

 1

		TO Method	CFD Solver	Sensitivity analysis	Optimization Algorithm	Application
	[82]	Density	FEM	Adjoint Variable Method	MMA	Cooling fins
	[83]	Density	n/a.	Discrete adjoint	SNOPT	Mono-fluid heat exchanger
	[84]	Density	FEM	Discrete adjoint	MMA	Micro-cooling device
	[50]	Density	FEM	n/a.	MMA	Mono-fluid heat exchanger
	[85]	Density	FEM	n/a.	MMA	Hierarchical microchannel
	[86]	Density	FEM	n/a.	MMA	Jet impingement target surfaces
						Multipass branching
	[87]	Density	FEM	n/a.	MMA	microchannel heat sink for electronics
						cooling
	[88]	Density	FEM	Adjoint Variable Method	MMA	Heat dissipating structure
	[89]	Density	FEM	n/a.	n/a.	Electro-fluid thermal compliant actuator
	[90]	Density	FEM	Continuous adjoint	SLP	Heat sink device
	[91]	Level set	FEM	Continuous adjoint	-	Mono-fluid heat exchanger
	[92]	Level set	FEM	-	Genetic algorithm	Mono fluid heat exchanger

1: Summary of applications of topology optimization on heat and mass transfer problems Table

  Une autre caractéristique importante de la fonction d'interpolation qui doit être discutée est le point stationnaire à ayant . En effet, l'algorithme d'optimisation complet qui cherche la valeur optimale de dépend de la dérivée locale au voisinage de la variable de conception sur la fonction d'interpolation. Ainsi la pente nulle pour représente une grande difficulté pour l'algorithme de franchir ce point pour faire une transition du fluide 1 au fluide 2 et vice versa. Ainsi la solution finale sera fortement dépendante de l'estimation initiale.

						est la suivante : pour le
	fluide 1	, pou le solide	et pour le fluide 2	. Ainsi pour	on aura un
	milieu poreux fluide 1-solide et pour	on aura un milieu poreux solide-fluide 2.
	Finalement la fonction d'interpolation est exprimée en (2.12) et les fonctions contraintes en
	(2.13) et (2.14).
	Pour assurer la convergence des variables de conception dans toutes les cellules vers 0, 0.5
	ou 1 et omettre toutes les valeurs intermédiaires de , la fonction d'interpolation avec laquelle le
	coefficient de perméabilité est calculé doit avoir une pente négative continue pour
	et	pour que la valeur minimale	soit strictement à	et	. Ainsi, cette
	fonction d'interpolation ne doit pas diminuer rapidement vers son minimum au voisinage de

, pour éviter d'avoir des cellules intermédiaires dans la solution finale ce qui nécessite une valeur comme le montre la Figure 2.4.

  control volume represented in Figure2.10 where its central node is P, its west, east, north and south faces are respectively w, e, n and s and its W, E, S and N are its neighbor nodes. It is then represented as follows:

																				23)
				Figure 2.10: Representation of a	control volume
	where subscripts		and		denotes respectively for			and	fields. Equation (2.22) is
	integrated over a  e e		w w	 	n n		s s		p P		u
	e		E	P		w		P	W		n		N	P		s		P	S	

  [START_REF]Plan Climat[END_REF]

	a b									a b
		a b	a b	a b	a b	a b	a b	a b	a b
	p									p
		w	e	s	n	e	w	n	s
	i j	i j	i j	i j	i j	i j	i j	i j	i j	i j
	a b	a b	a b
	w	w	w
	i j	i j	i j
	a b	a b	a b
	e	e	e
	i j	i j	i j
	ab
	s

Table 3 .

 3

				1: Continuity function coefficients at extremum values of and
	where , , and	are expressed as follows:
			11
	,	,	1,	1,	, 1	, 1	,	1,	1,	, 1	, 1
		1	2	2	2	2	2	1	1	1	1
			44
		11
					22

Table 4 .

 4

	1	1.3783	6.6031	36.37	163.63	127.26
	0.95	1.5667	7.1935	38.02	161.99	123.97
	0.94	1.747	7.629	39.86	160.13	120.27
	0.92	1.8735	7.8654	41.24	158.75	117.51
	0.88	2.2234	8.4958	43.76	156.23	112.47
	0.82	2.7179	8.9403	45.48	154.51	109.03
	0.7	3.5407	9.3732	47.55	152.44	104.89
	0.6	7.3929	10.429	53.18	146.81	93.63
	0.45	12.398	11.182	56.97	143.02	86.05
	0.3	18.633	11.744	59.53	140.46	80.93
	0.15	35.644	12.316	61.85	137.66	75.81
	0.05	41.327	12.402	62.48	136.45	73.98

5: Thermal and hydraulic performance of topology optimization results for various values of

(

Table 4 .

 4

	1	1.42	6.85	38.31	161.69	123.38
	0.94	1.58	7.77	42.26	157.73	115.47
	0.88	1.81	8.45	45.76	154.24	108.48
	0.8	2.12	8.94	50.56	149.38	98.82
	0.7	4.4	10.34	52.36	147.54	95.18
	0.55	8.64	11.67	58.69	141.38	82.69
	0.4	12.24	12.55	62.92	137.02	74.1
	0.3	15.67	13.04	65.65	134.19	68.54
	0.21	32.68	14.27	69.9	130.69	60.79
	0.15	33.88	14.35	72.41	127.04	54.63

7: Thermal and hydraulic performance of topology optimization results for various values of weighting factor (

  So it's clearly shown that for both cases the maximum heat transfer is expressed similarly. Hence heat exchanger effectiveness can be expressed as follows:

		q		c C T c o		c i T			C	min	, h i T		, c i T	(A2.1)
	For	we will have		and		:
		max q		 C T h h i ,		, h o T			C	min		, h i T		, c i T		(A2.2)

  Finally the total entransy dissipation number is the sum of equations (B2.3) and (B2.4) divided by the maximum entransy dissipation in a heat exchanger for non-dimensionless expression:

	P  G		  	m	P  	  	h	h o T	h o T T h i 		  	m	P  	  	c	c o T	c o T T c i 	(B2.4)
											T								T
											h i								c i
																			u
	w			W				P		W	WW					w
						G		 GG ,, h i c i TP Q T T   				(B2.5)
																		u
	e			P				E			P	W							e
																			u
	s			S				P			S	SS						s
																		u
	n			P				N		P	S							n
																			u
	w			P				W		P	E						w
					1
						3			2					for	u	0
	e			E				P			E	EE						e
					8

  3 are expressed in function of convective coefficients , diffusive coefficients and the source terms and defined respectively the coefficients are defined as follows:

						(					) (	)
	P	W	E	S	N	E	W		N	S
	W	W	W W
	S	S	S S
	E	E	 1	E		E
	N	N	 1	N		N
	 11  3 2 88 P W WW W W		3	W	2	P	E	 1	W		W
	 11  2 3 88 W P E E E			2	E	EE	3	P	 1	E		E
	 11  3 2 88 P S SS S S		3		S	2	P	N	 1	S		S
	1

Pour cela, trois méthodes différentes ont été suggérées dans ce chapitre, pour assurer la séparation des fluides :1. En introduisant une nouvelle fonction objectif qui pénalise la présence de deux fluides différents dans la même cellule, ou dans des cellules voisines. 2. En prenant l'effet de pénalisation comme une fonction de contrainte. 3. En modifiant le coefficient de perméabilité inverse (paragraphe 2.2.1) où le fluide 1 sera un obstacle pour le fluide 2 et vice versa, ce qui affectera directement la fonction objectif de dissipation de puissance totale du fluide .Afin de sélectionner la technique la plus appropriée, les trois méthodes ont été testées et comparées pour divers cas d'études.

solution presented in Figure 3.14 after 365 iterations. In comparison with the first two methods, fluid separation using constraint function is extremely costly in time, especially for big size mesh grids. Finally it is important to mention that the first variation in field occurred after 157 iterations, which means the method begun to separate fluids after the shape of pipes is approximately completely known. Thus, for configuration similar to double diffuser (paragraph 3.2.2.2) in which a conflict exist between fluid separation and absolute minimal fluid power dissipation in the domain, current method cannot prevent fusion of both pipes. The topology optimization problem described above was studied in mono-fluid topology optimization in [START_REF] Marck | Topology Optimization of Heat and Mass Transfer Problems: Laminar Flow[END_REF], [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF], [START_REF] Guest | Topology optimization of creeping fluid flows using a Darcy-Stokes finite element[END_REF], [START_REF] Sigmund | Topology optimization for multiphysics problems: A future FEMLAB application?[END_REF]. Mono-fluid results in literature showed that for an aspect ratio equal to 1, the optimal structure is two parallel straight pipes, similarly to the results in paragraph 3.2.2.1. However, for an aspect ratio AR=1.5, to get the lower fluid power dissipation, it is better to join the two pipes in a single wider pipe in the middle of the domain and to disjoin . The pressure drop due to bending section in optimal structure of , which has a , is . The geometrical parameters of all structures of between 1 and 0. For , reducing solid thickness between hot and cold streams for longest distance possible inside the optimization domain, doesn't satisfy anymore the required heat transfer rate. Hence, hot and cold pipes become having wavy shapes in order to increase heat transfer rate, and residence time of both fluids inside the domain by increasing the distance each fluid particle has to cross to leave the domain. It could be clearly seen from Figure 4.21, that when decreases, the length of cold and hot pipes increases and their diameter decreases, both lead to higher heat transfer rate and higher pressure drop. However, Table 4.1 shows that the increasing rate of for structures of is approximately similar for structures of , whereas the increasing rate of is much higher for structures of than for those of . This means when decreases, a high increasing in pressure drop is required in order to increase heat transfer rate by a small quantity. This could be also seen from Pareto frontier sketched in

Finally for extremely lower values of (

), the optimal structure is characterized by a chaotic behavior in optimizing the topology of cold and hot pipes, as seen in Similarly to parallel flow, the structures generated for the maximization of heat transfer for different values of are classified under two main categories:  For , the maximization of heat transfer is attained by bringing hot and cold pipes close to each other. The geometrical feature that characterizes these Minimization of power dissipation is then coupled with the maximization of heat transfer rate using multi-objective optimization technique, and continuity objective function which allowed designing a fluid to fluid heat exchanger device using the topology optimization technique.

However despite the success of the method proposed to deal with the problem posed at the beginning of the thesis, and to establish the basis of a design methodology of heat exchangers specifically, and various engineering devices that includes heat and mass transfer generally, the method still present some considerable limitations expressed as follows:

 The extremely high number of variables, which is double of the number of design variables for the same configuration in mono-fluid domain, limits the application of the method on relatively small domains, or force to use coarse meshes for large domains which decreases the accuracy of the final optimal structure, and decreases also the accuracy of the solution of CFD simulations at each iterations which dramatically affect the convergence process and the optimal result.  The convergence of the method to a reasonable solution is highly dependent on the numerical parameters of the various methods used in the optimization algorithm. The most important are the penalization parameters used in the interpolation function, the filter radius of sensibility and design fields, the step size in the method of moving asymptotes, the weight of continuity function in case of fluid separation using the method of objective function, and the value of in case of fluid separation using the modified inverse permeability coefficient. Usually those parameters vary dynamically throughout the optimization to ensure a smooth convergence and to avoid local optimum. However the choice of those parameters and their variation is so delicate which make the method extremely dependent on using the right set of parameters values throughout the optimization process.

Perspectives

Tridimensional domain

The first essential upgrade that should be performed on the method presented in this thesis is to consider the organization of materials in tridimensional domain, which will allow the method therefore to be considered as a new conceptual design methodology that leads to heat exchangers used in real life engineering applications. In fact the results of the bi-fluid topology optimization method applied on bi-dimensional domain could be considered as particular plan cut from the real 3D solution, while the optimization in 3D domain will lead to more complicated and realistic fluid paths that could not be predicted in 2D simulations. However 3D configurations increases the number of design variables, which will make the optimization problem more difficult due to the extremely high number of variables, which is already considered high in 2D simulations. This may require the use of more advanced numerical optimization algorithms to deal with the increment in number of variables, otherwise the 3D topology optimization method will be limited to extremely small domains or coarse mesh grid Appendices Annex C1: QUICK scheme Leonard [126] developed in 1979 a Quadratic Upstream Interpolation for Convective Kinetics (QUICK) scheme that computes cell face values using three neighbors nodes. The interpolation function is a quadratic function involving two neighbor nodes from each side of the face and a third node on the upstream side. For example suppose in Figure A where and are the bracketing nodes at each side of the face and is the upstream node that depends on the flow direction. Hence the QUICK scheme can be summarized as follows: