Somewhere, something incredible is waiting to be known.

Carl Sagan

Que pour examiner la vérité il est besoin, une fois dans sa vie, de mettre toutes choses en doute autant qu'il se peut.

René Descartes iii

Résumé en français

Les capteurs radio fréquence (RF) embarqués pour des applications aéronautiques ou spatiales ou pour les systèmes de communication doivent fonctionner à des puissances toujours plus élevées. Dans le domaine de la défense, ces niveaux de puissance sont nécessaires à la détection de cibles furtives ou afin de faire face aux brouilleurs. En ce qui concerne les systèmes de communication, de fortes puissances permettent d'assurer une bande passante suffisamment large pour faire face aux perturbations extérieures. Dans le même temps, ces capteurs doivent être de mieux en mieux intégrés afin répondre aux exigences de réduction de taille et de masse nécessaires dans le cas des drones ou des satellites. Seules des innovations en terme de management thermique des modules de puissance pourront permettre de progresser sur ces deux axes antagonistes.

Récemment, un nouveau type de semi-conducteur a vu le jour : le transistor GaN (Nitrure de Gallium). Le GaN est un matériau attractif pour les applications d'électronique de puissance à haute fréquence du fait de sa large bande interdite, de la grande vitesse de ses électrons et de sa forte conductivité thermique. Les performances du GaN sont, de fait, meilleures que celles des semi-conducteurs en silicium dont la faible bande interdite réduit la température de jonction ce qui limite la puissance du composant. En effet, dans de telles conditions, l'énergie thermique est suffisamment forte pour que les porteurs de charges passent la barrière d'énergie. Malgré cela, les puces GaN possèdent un inconvénient : afin d'éviter toute inductance parasite et de profiter au maximum de leur capacité de commutation à haute fréquence, leur packaging est conçu de telle sorte qu'elles puissent être montées sur le PCB (Printed Board Circuit) afin d'offrir la meilleure connexion possible. Pour cela, la croissance des puces doit se faire sur un substrat en epoxy à faible conductivité thermique, ce qui nécessite une solution de management thermique adapté.

Afin de faire face à la forte puissance et à l'augmentation rapide de la température de la puce, des solutions de management thermique efficaces sont nécessaires. Cellesci doivent pouvoir s'intégrer dans des espaces particulièrement réduits ce qui limite les possibilités et exclut de fait les solutions de refroidissement actif qui requièrent une circulation de fluide particulièrement encombrante. Les solutions actuelles, principalement v RÉSUMÉ EN FRANÇAIS constituées d'un dissipateur thermique en alliage métallique (aluminium ou cuivre principalement), ne suffisent plus a supporter les niveaux de puissance des nouveaux composants si leur masse et leur volume doivent être réduits.

De nouvelles solutions combinant forte capacité thermique, forte conductivité thermique et densités plus faibles doivent être développées. Ces dernières années un type de composés se démarque : les matériaux à changement de phase (PCM) qui sont des alcanes possédant une forte enthalpie de fusion permettant d'absorber de grandes quantités d'énergie au moment de la transition de phase. Différents PCM sont disponibles avec des températures de fusion couvrant une large gamme ce qui permet de s'adapter à une multitude de situations. Malgré ces nombreux avantages, les PCM possèdent une faible conductivité thermique qui limitent grandement la diffusion de chaleur, ralentissant transition de phase et absorption de chaleur.

Ce dernier point a mené au développement de matériaux composites, combinant des PCM à des structures possédant de fortes conductivités thermiques. Un vaste choix de renforts est disponible, ceux-ci se caractérisent par leur faible densité, cela va des poudre de carbone aux mousses métalliques, permettant d'améliorer la conductivité thermique tout en maintenant une forte proportion de PCM dans le composite, maximisant ainsi les capacités de stockage de chaleur. Bien que l'ajout de renforts soit absolument nécessaire, sa présence rend les simulations plus complexes. De fait, en plus de la simulation de la fusion, phénomène hautement non-linéaire, la conduction au sein du renfort thermique ainsi que l'échange thermique entre ce dernier et le PCM doivent être pris en compte. Une représentation complète du renfort demanderait une puissance de calcul trop élevée et des méthodes, dites d'homogénéisation, sont utilisées.

Cette étude doit permettre de déterminer les méthodes et hypothèses adaptées pour la simulation de dispositifs de management thermique pour l'électronique de puissance utilisant des PCM. Une fois validées expérimentalement, celles-ci doivent permettre d'évaluer les performances de nouveaux renforts thermiques fabriqués par impression 3D et dont la topologie sera optimisée. Ces structures doivent permettre d'obtenir de meilleures performances que les renforts actuels tout en permettant la fabrication des dispositifs d'un seul bloc, évitant de possibles résistances thermiques résiduelles.

Cette thèse CIFRE a été réalisée dans le cadre d'un projet de collaboration entre l'Institut d'Electronique, de Microélectronique et de Nanotechnologie (IEMN), de Thales LAS et de l'UMI CINTRA de Singapour (CNRS, Thales, Nanyang Technological University).

Le premier chapitre de cette thèse présente de façon détaillée les matériaux à changement de phase, leur composition ainsi que des détails sur leur première utilisation répertoriée (due à Maria Telkes et Eleonara Raymond) pour la Dover House. La suite du chapitre vi RÉSUMÉ EN FRANÇAIS détaille les différentes méthodes permettant de modéliser la fusion : des méthodes variable grid initiées par Stefan aux méthodes variable grid développées par Voller et al. dans les années 1990. Les premières permettent un suivi précis du front de fusion mais demandent une puissance de calcul supérieure, tandis que les secondes permettent une connaissance globale de la distribution de température de façon plus rapide. La suite du chapitre se concentre sur les renforts thermiques : le choix le plus adapté selon l'utilisation, ou encore la façon de modéliser leur comportement thermique mais également la façon de modéliser au mieux la relation de ceux-ci avec le PCM notamment en se concentrant sur les effets de Thermal Lagging.

Le deuxième chapitre présente le développement et la modélisation de nouveaux renforts architecturés fabriqués par impression 3D. En partant du constat que certaines propriétés des mousses métalliques (renfort utilisé majoritairement), notamment celles liées aux aspects aléatoires des mousses, avaient tendance à dégrader leurs performances, nous avons développé un type de renforts basé sur une structure périodique (lattice métallique). Dans une première partie, nous détaillons le développement d'un modèle de conductivité thermique effective (ETC) qui permet, de façon théorique, de montrer que le changement de topologie du renfort augmente de façon significative la conductivité effective. En comparant ce modèle à celui existant pour les mousses, une augmentation allant jusqu'à 70% est observée. La deuxième partie se concentre sur les effets de bord qui peuvent apparaître lorsque l'homogénéisation est réalisée avec un nombre d'unité de répétition trop faible (i.e. < 10). Ceux-ci sont dus à la constriction thermique qui apparaît lorsqu'un flux de chaleur se contracte du fait d'un changement brutal de géométrie.

C'est ce qui se passe lorsque l'on considère ce renfort dans son environnement global au contact avec le reste du dispositif. L'étude des travaux de Gladwell et al. a permis de déterminer l'influence de la constriction sur la conductivité réelle de ce renfort et de l'intégrer au modèle. Suite à cela, une série de mesures a permis de valider de façon expérimentale les modèles de conductivité.

En utilisant le modèle décrit dans le deuxième chapitre, nous développons un dispositif utilisant un renfort architecturé combiné à un PCM s'appliquant à un sujet d'étude concret dans le troisième chapitre. Pour cela, une série d'études paramétriques est réalisée afin de déterminer les paramètres les plus adaptés en termes de géométrie globale du dispositif, du choix du PCM et de topologie du renfort. Il apparaît ici qu'il existe un optimum de proportion de PCM dans le dispositif. Cette proportion permet que suffisamment de matériau conducteur (i.e. aluminium) soit présent dans le dispositif afin de dissiper efficacement et rapidement la chaleur mais également qu'une quantité suffisante de PCM soit présente pour absorber au mieux celle-ci. Afin de confirmer les résultats des simulations réalisées lors de l'étude paramétrique, une série de dispositifs est I would like to thank first Thomas Merlet (Thales LAS), Philippe Pernod (Centrale Lille/IEMN) and Philippe Coquet (Université de Lille/UMI CINTRA) for settling a global framework of collaboration between Thales LAS, the AIMAN-FILMS group of IEMN and UMI CINTRA. I would like to thank them, together with Olivier Bou Matar (Centrale Lille/IEMN), Jerome Foncin (Thales LAS), Dunlin Tan (Thales RT Singapore) and Edwin Teo (NTU) for establishing this thesis project within the latter framework. Then, thank to both Thales LAS and ANRT (Agence Nationale pour la Recherche et la Technologie) for their financial support.

For the past three years, I have been lucky to work alongside truly smart and dedicated people. I would like to thank Philippe Pernod who, despite constantly being on the rush, has found ways to help me in any way he could. Then, I have to thank Olivier Bou Matar who always found time to be there any time I needed a hand, thanks for the help and encouragements. Many thanks to Jerome Foncin for all the tips and precious knowledge you have conveyed along the way. Finally, I would like to thank Thomas Merlet for being such a caring and encouraging boss.

Then I would like to thank all the people I met at IEMN and who made that PhD such a great experience. Many thanks to Pierre, Manon, Aurélien, Thomas, Munique, Nicolas, Yannick, Abdelkrim, Stefano, Ghizlaine, Mathias, Laureline, Anna, Gong and I must forget some. Special thanks to Djamila for taking so much time for me.

I have to thank all people from Thales that have helped in any way. Many thanks to Benjamin, Simon, Paul, Alexis, Gauthier, Wafa, Arnaud, Julien, Anne-Charlotte, Guillaume, Robin and so many others. Time at Thales LAS would not have been the same without you guys. I have to thank Martin Raynaud from Thales Alenia Space too, for its always very interesting insights on the project.

I cannot forget all the people outside work, family and friends. Thank you to my parents for always have been so supportive in whatever I have embarked in. Many thanks to all my friends, the long term and the ones I met on the way, you all have always been there when I needed it the most.

Finally, special thanks to Julia for being such a wonderful person. Your help, at the very end of this work has been essential and I cannot thank you enough for that. Lord Kelvin (Tetrakaidecahedron) [START_REF] Thomson | On the Division of Space with Minimum Partitional Area[END_REF], (b) Weaire and Phelan [START_REF] Weaire | A counter-example to kelvin's conjecture on minimal surfaces[END_REF] and (c)

Gabbrielli [START_REF] Gabbrielli | A new counter-example to Kelvin's conjecture on minimal surfaces[END_REF]

General Introduction

Embedded Radio Frequency (RF) sensors for aircraft or space platforms, as for communication systems are more and more requesting powerful operating modes. For defense, high transmit power is for instance needed to cope with stealth targets and jammers and for communication systems, high transmitting power is key to ensure large bandwidth channels overcoming disturbances from adjacent radiations. Meanwhile, all these sensors also need a better integration to cope with size and weight reduction, sometimes mandatory as for Unmanned Aerial Vehicle (UAV) or satellites. These two axes of progress, inherently antagonists, shall only be addressed by breakthroughs as for the thermal management of the powerful modules.

In recent years, a new type of semi-conductors has emerged : GaN (Gallium Nitride) transistors. Actually, GaN is an attractive material for high power and high frequency devices due to its wide band gap, high electron velocity and high thermal conductivity.

It outperforms Si-based semiconductors that are limited by their relatively narrow band gap which lowers the maximum junction temperature hence maximum power. In fact, for such conditions, thermal energy is high enough for the charge carriers to pass the energy barrier. However, to avoid parasitic inductance and take full advantage of their high frequency switching capabilities, GaN chips are usually packaged so that they can be directly mounted on a Printed Circuit Board (PCB) to offer a maximum interconnection possibilities. To do so, the GaN chip must be grown on an epoxy substrate having a low thermal conductivity. This configuration necessitates thermal management solutions to prevent the GaN chip from overheating.

To cope with the high power and high temperature rise in the chip, efficient thermal management solutions are required. These solutions must be integrated in tight space environments which limits the possibilities and usually exclude the use of active thermal management solutions requiring fluid circulation, for example. Current solutions, in embedded RF sensors, involve heat sinks or heat spreaders mostly made of metals. Although they have been proven reliable, such solutions can no longer cope with the higher power density and mass and size reduction requirement.

New solutions with higher thermal capacity, high thermal conductivity and lower den-

GENERAL INTRODUCTION

Figure 1: PCM melting in an idealized metal foam filler.

sity needs to be found. Over the past decades a solution has emerged : Phase Change Material (PCM) are alkane-based materials possessing a high latent heat of fusion allowing the absorption of large heat quantity during phase transition. They are available with melting temperatures spanning over a large range that can be chosen to meet a variety of situations. Despite their advantages, PCM suffer from their low thermal conductivity which limit heat spreading hence phase transition and heat absorption.

For this reason, composite materials combining PCM and highly conductive structures are being developed. A wide variety of fillers are available in the literature ranging from carbon-based powders to metal foams which allows thermal conductivity enhancement while maintaining a high proportion of PCM in the composite, maximizing heat storage.

While the addition of the filler enhances the PCM, its intricate geometry makes it more difficult to simulate it. Actually, in addition to the fusion, a highly non linear phenomena, heat conduction in the filler as well as heat exchange between the PCM and the filler itself must be accounted for (see figure 1). Representing the entirety of the filler would be time consuming and so-called homogenization techniques are preferred.

The present work aims at determining the right methods, as well as the right assumptions to perform accurate and efficient simulations describing the behavior of PCM based device in the context of power electronic thermal management. Once set and validated experimentally, simulations must help assess the performances of new types of filler based on deterministic topology and additively manufactured. Such structures are thought to outperform traditional fillers as well as facilitating integration as devices are manufactured as a single block avoiding any thermal resistance that could limit their performances.

The present manuscript is composed of four chapters detailed below :

• Chapter 1 details the use of PCM as found in the literature. First, a definition of the PCM is proposed. Then, the chapter concentrates on PCM modeling with emphasis on fusion modeling and interactions with the filler. Finally, different ways to enhance PCM are proposed as well as way to assess each solution's performances.

GENERAL INTRODUCTION

• Chapter 2 introduces a new type of PCM filler based on the use of lattice structure which topology should, for a given conductive material volume ratio, outperform traditional filler performances. To prove it an Effective Thermal Conductivity (ETC) model is developed and compared to the existing ones for traditional (i.e. foam) filler. While it was first done theoretically, considering a single isolated cell, the lattice filler is then considered in its actual environment which reveals side effect due to constriction phenomena. These phenomena could, in some case, drastically reduce the ETC of the resulting filler which is detrimental to the performance of the PCM based device. In a last part, the previously described models are experimentally validated.

• To better assess the performance of the filler developed in chapter 2, the ETC model is implemented in an actual use-case in chapter 3. Using the different modeling techniques from the literature and the present study, a device aiming at limiting the temperature rise of four 100 W dies is optimized in terms of package and PCM enclosure dimensions, filler properties and PCM choice, and compared to existing thermal management solutions. This theoretical study is then confirmed experimentally to validate the simulation assumptions.

• In chapter 4, a more complex and theoretically more efficient of PCM enhancement technique is assessed : Cascaded storage. Rather than using an homogeneous filler, the ratio of conductive material is modified depending on location in the device, usually more conductive materials close to the heat source provide better results.

To provide the optimum material distribution, a combination of a space colonization and a genetic algorithm is used. This allows the creation of a so-called heat spreader added in the PCM enclosure which shape is optimized to minimize die temperature. In the meantime, the porosity of the lattice filler around is adapted to maintain similar quantities of PCM maximizing heat storage and storage rate.

Introduction

As seen in introduction, as integration and power density of electronic components increase, thermal management become a vital part of any new product development. Solutions involving heat spreading or heat exchange have been widely used but lack of space or inability to efficiently exchange heat because of the electronics surroundings might require additional solutions. The use of Phase Change Material (PCM) has become a major topic of research in the last decades due to their high latent of fusion. The high latent heat allows for more energy storing in a smaller mass or volume compared to sensible heat based devices. In addition, using latent heat allows for energy storing at a constant temperature corresponding to the melting point of the PCM [1] which is of great importance while considering thermal management issues.

Despite their great latent heat of fusion PCM often suffer from their low thermal conductivity that usually limits the charging or discharging time. For better integration, especially when dealing with high power density applications, thermal conductivity enhancement is required. In that purpose a large amount of different fillers has been developed and used; one can cite : high-conductivity dispersed filler [2; 3], carbon-based [4][5][6], or metal [7; 8] foams or architected metal structures [9]. Those solutions have seen thermal conductivity enhancements ranging from 2 to about 15 (enhancement normalized to the ratio of filler in the composite [5]). Although improvements have been made in terms of thermal conductivity enhancements, such composite might not be as efficient as needed for some cases and additional solutions may be required. In that sense, several solutions involving a variation of material properties in the entire composite have been developed.

The variation of thermal properties in the composite allows for better heat spreading close to the heat sources while storing is performed away for them. This technique, usually referred to as Cascaded storage [10] are relatively diverse as they involve : several PCM of different melting point and latent heat in cascade [10], similar PCM infused in fillers of varying density and thermal conductivity [11][12][13], combination of finned heat spreader and enhanced PCM [14] or PCM placed in topologically optimized heat spreader [15; 16].

Assessing the efficiency of a given PCM based solution can be done experimentally but this could turn out to be costly, especially if the solution needs to be optimized for a given issue. In that sense, efficient simulation techniques are required to simulate fusion which is a highly non-linear phenomenon. For that, two different strategies have been adopted. The first one (also called Stefan problem [17]) involves two distinct materials with a moving boundary that is accounted for through mesh rearrangement at each time step, which is computationally costly. On the other hand a so-called fixed grid problem is defined [18], involving only one material with a varying specific heat capacity that accounts for the additional absorbed energy around the melting point. The latter does not CHAPTER 1. PHASE CHANGE MATERIALS require mesh rearrangement which save computation time. Going further in the assessment of a PCM based device requires to account for the heat exchange between the PCM and the conductive matrix in which it is embedded. To do so, two strategies are found in the litterature : Local Thermal Equilibrium (LTE) or Local Thermal Non-Equilibrium (LTNE) [19][20][21]. The former considers that continuity of temperature is verified at the interface between the PCM and the conductive matrix while the former considers a temperature gap and a convective-like exchange coefficient that is to be determined for each situation.

In the first paragraph of this chapter, the use of PCM as energy storing or thermal management solution is reviewed in terms of possible application and of PCM choice depending on the required melting temperature and latent heat of fusion. In the second paragraph, the issues of Fusion Modeling and Homogenization Theory are discussed. Finally, the third paragraph of the chapter concentrates on PCM enhancement techniques and particular attention is given to the modeling of foam as PCM enhancer and to Cascaded Storage.

Phase Change Materials Definition and Use

What is a PCM?

Although any material can experience phase change, whether it is a solid/liquid, liquid/gas, solid/gas or a crystalline structure transition, only some specific materials could be used as PCM for energy storage or thermal management solutions. To be considered a viable PCM, a given material must respect the following conditions [1; 22]:

• a melting point in the application range,

• a high latent heat of transition per unit mass to minimize the required mass of PCM,

• high specific heat that provides addition storing effect through sensible heat,

• high thermal conductivity to minimize charging and discharging time,

• limited thermal expansion, especially at the melting point, so that simple container can be used,

• high nucleation rate to avoid or limit subcooling during freezing,

• high rate of crystal growth to meet demands of heat recovery from the storage system,

• complete reversible freezing/melting cycle properties, • possess chemical stability, no chemical decomposition and corrosion resistance to construction materials,

• contain non-poisonous, non-flammable and non-explosive elements/compounds,

• available in large quantities at low cost.

Only some materials respect at least some of the above mentioned characteristics and can be divided into three categories (see figure 1.1) [22; 23] : Organic compounds, Inorganic compounds and Eutectics.

Organic compounds

Organic PCM can undergo phase transition over a large number of cycles as they are not subjected to phase segregation or degradation of their latent heat of phase transition.

Other advantages of such PCM include little or no super-cooling and low corrosiveness.

Organic compounds can be subdivided into two sub-categories: Paraffin and Non-Paraffin organic PCM.

• paraffins usually consist of a mixture of n-alkanes (i.e. CH 3 (CH 2 ) n CH 3 ). The crystallisation of the CH 3 part releases a large amount of energy that usually increases with the length of the chain. Paraffin are reliable, safe, cheap, non-corrosive and available in large quantities. Their melting point have been reported to range from -37 to 111 • C1 . In addition, they are usually stable below 500

Inorganic compounds

As for organic compounds, inorganic ones have the advantage of not degrading with cycling and not super-cooling. They can be divided into two subcategories as follow : Salt hydrates and Metallics.

• Salt Hydrates of general chemical formula AB • nH 2 O are alloys of an inorganic salt AB and of water added through crystallization.During phase transition, dehydration occurs leading to the formation of a less hydrated compound through the reaction:

AB.nH 2 O → AB.mH 2 O + (n -m)H 2 0, (1.1)
or to an anhydrous form:

AB.nH 2 O → AB + nH 2 O. (1.2)
Depending on the predominant transition reaction and anhydrous salt solubility in water, different behaviors are expected. As a matter of fact, congruent melting (i.e. the product of the reaction remains in the same phase) when the salt is soluble in water at melting temperature, while incongruent melting occurs when the salt is only partially soluble in water. A third way called semi-congruent melting occurs when solid and liquid phases in equilibrium have different compositions which usually happens when equation 1.1 is predominant [10]. Salt Hydrates have been extensively studied and the following general properties have been reported : higher thermal conductivity than Organic compounds (around 0.5

W•K -1 •m -1 , i.e.
double organic compounds thermal conductivity), higher heat of fusion than organic compounds (see figure 1.2) but they are usually corrosive and release water which, depending on experimental conditions, could turn out to be a threat to the device integrity. • Metallic PCM are composed of low melting point metals and some of their alloys such as: Aluminum, Magnesium, Zinc or Copper. They differ from other PCM by their high thermal conductivity (usually over 100 W•K -1 •m -1 ) and high density. While their high density might be seen as a drawback in applications in which mass reduction is a requirement, they could turn out to be the optimum solutions when volume reduction is sought. Metallics are particularly well adapted to high temperature and high power applications due to their high thermal conductivities and melting points.

Eutectics

Eutectics are a combination of two or more compounds that melts at constant temperature and composition. They always undergo congruent melting allowing highly reproducible melting/freezing cycles. Their overall composition are particularly diverse hence are their melting point and latent heat of fusion (see figure 1.2). Eutectics can be divided into three sub-categories depending on their compositions: organic-organic, inorganicorganic and inorganic-inorganic.

CHAPTER 1. PHASE CHANGE MATERIALS PCM Type T m ( • C) ∆H (J•g -1 ) k (W•K -1 •m -1 ) ρ (kg•m -3 )
• Organic-organic eutectics are combinations of two or more fatty acids [24; 25]. Compared to isolated fatty acids they have a lower melting point that the ones composing them which is more appropriate for, for example, heating/freezing of buildings.

Another advantage of binary eutectics made from fatty acids is their high surface tension that facilitates their impregnation in a porous thermally conductive matrix [24].

• Inorganic-organic and inorganic-inorganic eutecticts are similar compounds as they are the combination of a metallic or organic compound and of one or more halogen salts (Fluioride, Chloride, Bromide, Ionide...). Their melting points as well as enthalpy of phase transition are highly diverse [26] (see figure 1.2) but usually superior than other types of PCM.

In table 1.1, the characteristics of the different PCM detailed in the previous paragraphs are reported to better compare and adapt them to the desired application.

History of PCM as Energy Storage or Thermal Management Solutions

As described before, the main interest of PCM lies in their high latent heat of fusion which allows the storage of a large quantity of energy at an almost constant temperature. Added to a more complex system, this has a tendency to increase its Thermal Inertia. Thermal Inertia is defined as the "property of a material that expresses the degree of slowness with which its temperature reaches that of the environment" or as the "capacity of a material to store heat and to delay its transmission" [27]. In that sense, PCM main applications include, but are not limited to [28]: • Thermal damping: the PCM is used to limit the variation of temperature of a system placed in an oscillatory temperature environment. Such examples are found in the aerospace industry to limit the temperature variations due to sun exposition, or in the building industry to cope with the day/night cycle and ensure thermal comfort (see figure 1.3a).

• Inhibiting Thermal Excursion: the PCM is used to limit or delay temperature increase while the system is exposed to a brutal change in the environment temperature for a limited period of time. This is found in high power applications like batteries or more generally high power electronics (see figure 1.3b).

• Maintaining Constant Temperature: PCM is used here to maintain the system at a constant temperature (i.e. melting temperature) as long as it has not undergo a complete phase change. This only works for low power and relatively short term applications like food, beverages or drug transportation, heat and cold therapy for muscle and joint pain... (see figure 1.3c).

Use of solar energy pioneer work: the Dover House

The first ever mention of PCM as thermal energy storage solution is due to Maria Telkes (Pioneer of the Thermoelectricity [30; 31]) and Eleanor Raymond (Architect) [29; 32] who designed the first solar heated house called the Dover House in 1949 (see figure 1.4).

Telkes proposes here a first solution to avoid the use of fossil fuel for house heating and rely on solar energy. C). Heat from the collector is conducted to the bins through ducts and fans are used to ensure hot air circulation to the room to be heated as shown on figure 1.5. The system developed by Telkes includes 21 tones of PCM that can stores 4 × 10 6 Btu (equivalent to 4 × 10 9 J) which represents ten days of average heating requirements. This means that once fully charged, this system can heat the house for ten consecutive days. It was determined that this system was able to convert 41% of the solar radiation received each day into heat with peaks at 60% during sunny days. During summer the exact same system could be used to lower the temperature by insulating the collector from the outside. In this configuration, the PCM absorb heat from the room during daytime and the cool air at night helps release the heat in the atmosphere [29]. Although this project looked promising and received a lot of attention at the time it was realised, after the third winter, salt hydrates showed segregation and leaked off the bin due to corrosion [33] making the house inhabitable without fossil fuel heating.

Evolution and increasing interest for PCM

Due to the inherent difficulties linked to the actual application of the PCM, they did not receive much attention after Telkes and Raymond work due to the increasing offer of nuclear power and cheap oil from the Middle East [33]. This was true until the late 1970s and early 1980s when the world was hit by the first energy crisis. At the time, a renewed interest for the use of PCM in the building industry is noted [1] especially for the storage of solar energy. After the crisis, interest dropped again until it recently appeared as a so- lution to solve thermal management issues in a increased number of fields. For example it is thought to help solving the problem of fluctuation of offer and demand on the power grid as it could store the excess of power produced off-peak hours and release it during peak hours [34]. This is especially important for the production of electricity through renewable energy such as wind or solar [35] as the supply is extremely fluctuating. A second example is related to the thermal management of batteries which production has grown exponentially during the past decades due to the rise of electric or hybrid vehicles. Thermal management of battery pack is crucial as operation out of the optimal range defined by the manufacturer often limits the charging capacity and shortens the battery life span [36; 37]. Lazrak et al. [36] report a 30% life span reduction for batteries functioning only 5 • C over the temperature range defined by the manufacturer. Finally, PCM are thought to be a good fit to produce compact and light heat spreader for the thermal management of power electronics through smart thermal conductivity enhancement and packaging [14; 38; 39]. Such devices must cope with the increasing power density encountered in such situation (up to 100 W for commercial applications and up to 1000 W for military applications) combined with harsh environments (vibration, shock, high temperature...)

and lack of space.

The development of PCM as viable energy storage system, or Latent Heat Thermal Energy Storage System (LHTESS) as they are usually called, is only possible if the issues of fast charging/discharging, optimum average temperature, compactness, light weight and low cost are solved. To ensure they all are, the understanding of a number of physical phenomena that are still active research topics today are absolutely necessary:

CHAPTER 1. PHASE CHANGE MATERIALS
• Phase-change problem formulation to account for the moving boundary between solid and liquid [18; 40]. Most of the solutions involved the use of the so-called Enthalpy Method that describes the melting PCM as a unique material with a varying specific heat capacity avoiding the need to actually track the moving interface between liquid and solid.

• Thermal conductivity enhancement [41] to solve one of the major issue related to PCM which highly limits their performances. The composite created by the infusion of such filler in the PCM could turn out to have a complex geometry which is difficult to accurately represent in simulation software without considerable computing resources. Hence, numerous researches concentrated on the best way to model such filler and determine the ETC [START_REF] Yovanovich | Spreading resistance of isoflux rectangles and strips on compound flux channels[END_REF][43][44][45].

• While the modeling of ETC is related to the steady-state behavior of the PCM composite, its dynamic behavior involves heat exchange between the filler and the PCM itself. Due to the combined effects of the complex interface topology and the great difference in thermal properties, continuity of the temperature at the interface is not guaranteed. This subject, referred to as homogenization, discuss the possibility of describing the thermal behavior of such composite with one (One-Temperature Model or LTE) or two (Two-Temperature Model or LTNE) sets of differential equations and the limits of each representation [19; 20].

Fusion Modeling and Homogenization theory 1.3.1 Fusion Modeling : Variable against Fixed Grid Methods

As mentioned in the previous paragraph, accurately describing and modeling fusion is key to assess the viability of a PCM based solution for energy storage or thermal management. The first model used to describe fusion was developed by Stefan [17] in 1891. In his model the heat equation is solved on a domain Ω separated into two sub-domains representing the two phases of the material : liquid denoted Ω 1 and solid denoted Ω 2 .

The position of the moving boundary between the two phases is called s and depends on the the time t and the temperature at the boundary is the melting temperature of the considered material denoted u m [46; 47].

In its first, most simple formulation, the Stefan problem is defined on a one-dimensional semi-infinite region (0 ≤ x < ∞) with a melting temperature u m . In the following paragraph, the equations of the model are defined.

CHAPTER 1. PHASE CHANGE MATERIALS

In the liquid region Ω 1 (0 ≤ x < s(t )), the heat equation is solved:

∂u ∂t = α l ∂ 2 u ∂x 2 , (1.3)
with α l the liquid thermal diffusivity.

In the solid region Ω 2 (s(t ) ≤ x < ∞), the temperature, for simplification purpose, is considered constant and equal to the melting temperature u m , hence:

u(x, t ) = u m . (1.4)
To describe the different boundary and initial conditions, ones must first define the so-called Stefan condition that defined the heat flux at the solid/liquid interface. To do so, the domain Ω is represented as a cylinder along the axis x having a cross sectional area S.

Then, one must express the energy E required for the boundary to move from s(t 0 ) to s(t 1 )

with t 0 and t 1 > t 0 two different times:

E = S(s(t 1 ) -s(t 0 ))ρ l L, (1.5)
with ρ l the liquid density and L the material latent heat of fusion. This energy can also be expressed as the combination of the two heat fluxes coming from the solid (φ s ) and the liquid (φ l ) at the interface forming a disk A:

E = t 1 t 0 A [φ l n + φ s (-n)]d Ad τ = t 1 t 0 S -k l ∂u ∂x Ω 1 + k s ∂u ∂x Ω 2 d τ, (1.6)
with n a unit vector in the x direction, k l and k s are respectively the liquid and solid thermal conductivities.

Substituting equation 1.6 into equation 1.5 and dividing both terms by (t 1t 0 ) leads to:

s(t 1 ) -s(t 0 ) t 1 -t 0 ρ l L = 1 t 1 -t 0 t 1 t 0 -k l ∂u ∂x Ω 1 + k s ∂u ∂x Ω 2 d τ. (1.7)
If one considers now the limit of the previous equation when t 1 -→ t 0 equation 1.7 becomes:

d s d t t 0 ρ l L = lim t 1 →t 0 1 t 1 -t 0 t 1 t 0 -k l ∂u ∂x Ω 1 + k s ∂u ∂x Ω 2 d τ, (1.8)
which can be, using the Mean-Value Theorem for integral [46], written as: This condition is applied to the specific problem described earlier and the two boundary conditions at x = 0 and x = s(t ) are defined, being respectively:

d s d t t ρ l L = -k l ∂u ∂x Ω 1 + k s ∂u ∂x Ω 2 , ( 1
u(0, t ) = u 0 , (1.10)
and

d s d t t ρ l L = -k l ∂u ∂x s(t ),t , (1.11)
while the initial condition is simply written:

u(x, t ) = u m . (1.12)
Analytical solutions exist for this problem and are expressed as :

           u(x, t ) = u 0 1 -1 erf(λ) erf x 2 α l t s(t ) = 2λ α l t λe λ 2 erf(λ) = S t ,L π , (1.13)
with erf the error function and S t ,L = C p,l (u 0 -u m )/L the Stefan number in the liquid (C p,l is the specific heat capacity of the liquid) and λ is a dimensionless temperature.

Figures 1.6 depicts the evolution of the temperature distribution along the axis x and the evolution of s(t ) at arbitrary moments in the simulation. Although this method is relatively simple and gives a first approximation of the modelling of the fusion of a material, it cannot represent realistic situation as the temperature in the solid remains constant.

In addition, analytical solutions of the Stefan problem only exist for a limited number of boundary and initial conditions [48]. In addition it is adapted to compounds melting at CHAPTER 1. PHASE CHANGE MATERIALS a given temperature but not on a temperature range as it is the case for most materials, especially for mixture or composites which represent a large part of all PCM mentioned in the previous paragraph hence numerical methods are needed. In that sens two strategies have been adopted which are referred to as variable or fixed grid methods [47].

Variable Grid Methods

Solving a physical problem numerically requires to divide time and space into discrete points at which the problem will actually be solved while at any other position the value of the variable that is determined will only be derived through interpolation. Difficulties arise when this is applied to a fusion problem : one has no way to determine if, at the next time step, a point is liquid or solid hence making the fusion front tracking relatively complicated. To solve this issue a first set of methods has been developed in the 1980s called the Variable Grid Methods. Using them, either the time [49][50][51] or the space [52][53][54] steps are modified at each time step so that the melting front match the position of a given point of the grid.

The interface-fitting grids [47] are methods using a uniform spatial grid but a varying time step that were developed to solve two-phase and one dimensional problem similar to the Stefan problem but adapted to a higher number of boundary and initial conditions.

These methods solve dimensionless form of the heat equation.

∂u ∂t = ∂ 2 u ∂x 2 , (1.14)
for 0 ≤ x(t ) and t > 0 with the following boundary and initial conditions:

                 ∂u ∂x = -1 (x = 0, t > 0), u = 0 (x ≥ s(t ), t > 0), d s d t = -∂u ∂x (x = s(t ), t > 0), s(0) = 0. (1.15)
To conform with the discretization of the problem a discretized form of equation 1.14 using a backward difference scheme is defined with the dimensionless temperature u i , j , given at a time t = j ∆t and a position x = i ∆x at the k th iteration of the time-step determination algorithm by:

u (k) i , j +1 -u i , j ∆t (k) j = u (k) i +1, j +1 -2u (k) i , j +1 + u (k) i -1, j +1 ∆x 2 . (1.16)
To determine the proper time step ∆t j at the j th instant, Douglas et al. [49] integrated equation 1.14 with respect to x between 0 and s(t ) so that the position of the melting front CHAPTER 1. PHASE CHANGE MATERIALS actually match the position of the j th point of the spatial grid, leading to:

t = s(t ) + s(t ) 0 ud x, (1.17)
which could be discretized as:

∆t (k+1) j = ( j + 1) + j i =1 u (k) i , j +1 ∆x -t j .
(1.18)

A higher value of k only increases the precision of the temperature distribution and position of the melting but obviously increase computational cost as well.

Using the value of the time step, equation 1.16 can be solved numerically. Other methods defining the proper variable time step were developed by Goodling et al. [51] and Gupta et al. [50] which showed successive improvements. Gupta et al. compared their results and those from Douglass [49] and Goodling et al. [51] with an integral methods which showed "surprisingly good agreement" [50] confirming the interest of such method. However Hu [47] pointed out that this method can only be applied to onedimensional problems, limiting its use to actual engineering issues.

The other type of Variable Grid Methods are based on variable space grid also known as dynamic grids [47]. Using these methods, the total number of points in the grid is kept constant but the position of those points is adjusted at each time step so that the boundary s(t ) lies on a particular grid point. The method was first developed by Murray et al. [52] and further explored by Heitz [53] and Crank et al. [54]. To account for the moving grid, equation 1.14 needs to be modified. In this case, the time derivative of any point i of the grid is expressed as:

∂u ∂t i = ∂u ∂x x,t ∂x ∂t i + ∂u ∂t x,t , (1.19)
where the moving rate of the grid is directly related to the moving boundary location s(t ):

∂x ∂t i = x s(t ) d s d t . (1.20)
Substituting equations 1.19 and 1.20 into equation 1.14 leads to:

∂u ∂t i = x s(t ) d s d t ∂u ∂x x,t + ∂ 2 u ∂x 2 x,t (1.21)
that is solved using a finite difference scheme by Heitz [53] and Murray [52]. To do so, the space is discretized into N total points: n points in the liquid separated by a spatial interval ∆x l = s/n and Nn in the solid separated by an interval ∆x s = (Es)/(Nn) (E being the total length of the material). In the liquid, if the material is undergoing fusion, spatial interval increases with time while the contrary is true in the solid. Although this CHAPTER 1. PHASE CHANGE MATERIALS method properly works for most moments of the simulation, if only one phase is present at any moment (i.e. s(t ) = 0 or s(t ) = E) equations 1.20 and 1.21 are not defined which requires the definition of initial assumed value s 0 and a starting temperature distribution and to stop the simulation before s(t ) = E, extrapolate the temperature distribution at that particular moment and then go on with the following of the simulation [52]. In that sens, additional methods involving Taylor expansion [52], Cubic Splines [54] lead to better results [47]. Although these methods showed great agreement with experimental results

for one-dimensional problem, their extension to two or three-dimensional problems requires precise tracking of the fusion front at several spaces of the domain demanding additional computational power.

Fixed Grid Methods

When a precise determination of the fusion front location is not required but focus is rather made on the temperature distribution evaluation for multi-dimensional problems, the use of a fixed grid method is more adapted [18; 40; 48; 55]. Unlike variable grid methods, in which fusion and the related heat exchange happen at the interface between two different phases, fixed grid methods consider only one phase for which fusion is accounted for through equations of state. Most fixed grid methods are based on the same principle arising from the so-called Enthalpy Method. Using this method the heat equation is reformulated in term of enthalpy H rather than temperature only:

k∇ 2 T = ρ ∂H ∂t , (1.22)
The enthalpy term H is expressed as the sum of two terms accounting for both sensible heat h and latent heat ∆H:

H = h + ∆H, (1.23)
with:

h(T) = T 0 C p d T, ∆H(T) = L f l (T). (1.24)
in which C p is the specific heat capacity of the compound, considered unchanged after phase transition and f l ∈ [0, 1] is the liquid fraction at a given point of the domain Ω.

The expression of the latter depends on the nature of the fusion or solidification (kinetics, subcooling, nucleation rate, segregation...). Used as so, this formulation is robust and relatively insensible to space and time steps. However, König et al. [56] pointed out the step-behavior of such method as the temperature evolution with time may show some CHAPTER 1. PHASE CHANGE MATERIALS plateau depending on the Stefan number and the thermal conductivity differences between the materials found in the simulation. This is due to the fact the nodal enthalpy, an average value, is used as a local thermodynamic enthalpy.

Equation 1.22 is non-linear, as it contains two related but unknown variables H and T.

It could be convenient in certain case to reformulate this governing equation in terms of a single unknown variable T [55] with a non-linear latent heat effects "isolated" in a source term coefficient. To do so, equation 1.22 is modified as so:

k∇ 2 T = ρ ∂T ∂t ∂H ∂T , (1.25)
The derivative of H according to T could be interpreted in different ways that may change the way Finite Element Analysis (FEA) algorithms solve the physical problem. The first interpretation uses equation 1.24 which leads to:

k∇ 2 T = ρC p ∂T ∂t + ρL ∂ f l ∂t . (1.26)
which is called the Latent Heat Source Term Method. In that situation, the heat equation is not written in terms of a single variable as f l appearing in the source term is also an unknown function of T. Updating the value of the f l field is actually key feature in obtaining a reliable numerical solution. The common representation of f l is reported on equation 1.27 although others, with less physical meaning can be found [56; 57]. The determination of the actual value at each point of Ω for each time step is done through iterative methods [55]:

f l =            0 T < T s , T-T s T l -T s T s ≤ T < T l , 1 T ≥ T l .
(1.27)

with T l and T s respectively the liquidus and solidus temperatures.

This method is probably one of the most used [14; 15; 39; 58-60] as it gives results as reliable as the enthalpy method with similar default. However, the convergence of the iterative method may not be ensured for isothermal transition, convection-dominated phase change or when high temperature gradient are encountered.

The enthalpy derivative in equation 1.25 could also be interpreted in a different way.

Actually, this term is the definition of the specific heat capacity from a thermodynamic point of view and is used to defined a non-linear parameter called Apparent c app or Effective Heat Capacity c eff which both vary with temperature to account for the phase change.

Although, those two terms refer to two slightly different methods, their mathematical formulation is identical as equation 1.25 becomes:

k∇ 2 T = ρc i ∂T ∂t .
(1.28)

CHAPTER 1. PHASE CHANGE MATERIALS

The only difference between those two methods is the definition of the heat capacity:

using the apparent heat capacity method, heat capacity is determined in terms of the nodal temperature while the effective heat capacity is calculated as an average value on a control volume of volume V which tends to increase precision, hence:

c eff = 1 V V c app d V. (1.29)
The main advantage of both methods is their easy implementation in one-phase heat transfer problem. However, these methods are troublesome if the temperature rise is high enough so that temperature of a given point rise from a temperature below the solidus to above the liquidus. In this case, the energy absorption due to the phase change might not be accounted for. These methods are also known for resulting in temperature oscillation which can be solved by time and space steps reduction which could be computationally costly.

Homogenization Theory: Local Thermal Equilibrium (LTE) or Non-Equilibrium (LTNE)

In the previous paragraph, ways to model PCM fusion have been listed and detailed which shows the potential of the enthalpy method when any material is considered by itself.

However, as said in part 1.2.2, PCM need to be infused in a conductive filler to enhance their thermal conductivity and overall performances as LHTESS. Although it has been done only on parts of such composite [14; 58; 61], representing the entirety of the filler might increase computational cost as most filler are foams, architected structures or powder (see part 1.4) having complex structures that obviously increase the number of points in the space grid. Simulating such composite this way is called Direct Numerical Simulation (DNS) and requires at least two sets of differential equations if all convective effect are neglected and if the apparent capacity method is used, to describe its thermal behavior :

ρ f C p, f (T f ) ∂T f ∂t = ∇ • (k f ∇T f ), (1.30)
in the liquid, and:

ρ s C p,s ∂T s ∂t = ∇ • (k s ∇T s ), (1.31)
in the solid. In equations 1.30 and 1.31, subscript f refers to the PCM (or fluid) properties while s refers to the filler (or solid) properties. To those equations must be added the interfacial thermal boundary conditions:

T f = T s , k f ∇T f • n = k s ∇T s • n.
(1.32)

CHAPTER 1. PHASE CHANGE MATERIALS at the PCM/filler interface.

To reduce computational cost, so-called Volume-averaged simulations in which the composite created is represented as bulk material having properties that depends on both filler and PCM properties have been developed. The most simple volume-averaged model is called the one-temperature model based on LTE assumption [58]. Using it, it is assumed that there is no difference between the fluid and the solid respective average temperatures simplifying equations1.30 and 1.31 as they can now be written as one unique equation, using an average temperature 〈T〉 written:

〈T〉 = 1 V V Td V, (1.33)
in which V is the total volume of the domain Ω in which the problem is solved. Using the average temperature, it leads to:

ρ eff C p,eff ∂〈T〉 ∂t = ∇ • (k eff ∇〈T〉), (1.34)
with the subscript eff referring to the effective properties depending on both fluid and solid properties and on the fluid volume ratio ε. While the expression of k eff will be discussed later and is key to obtain a precise result, the expression of ρ eff and C p,eff are defined by the following mixing laws :

ρ eff = ρ f ε + ρ s (1 -ε), (1.35)
and

C p,eff = ρ f C p, f ε + ρ s C p,s (1 -ε) ρ f ε + ρ s (1 -ε) . (1.36)
Although this definition is relatively simple and has been proven to be valid for a variety of cases [12; 14; 58; 61; 62], it was pointed out by Quintard et al. [19][20][21] that if the characteristic length of the network formed by the solid filler or if the thermal conductivities the filler and the PCM were too different, error may arise as the LTE assumption might not be valid anymore. This observation was originally made by Cattaneo and Vernotte [63][64][65] who expressed the impossibility of using the traditional heat equation in low pressure gases in which the instantaneous transmission of heat arising from Fourier law was not verified anymore. In such case, the average distance between two particles is higher than their mean free path which proscribe instantaneous heat transfer. Consequently, Cattaneo and Vernotte provide a correction to the heat equation and Fourier law introducing a second time derivative of the temperature and a time τ which is characteristic to the mean time between two particles collision. In that case the Fourier law becomes:

-k∇T = q + τ ∂q ∂t , (1.37)
CHAPTER 1. PHASE CHANGE MATERIALS with q the heat flux. Hence, the heat equation is modified as so:

α∇ 2 T = ∂T ∂t + τ ∂ 2 T ∂t 2 , (1.38)
with α the thermal diffusivity.

This equation has been used to solve a number of heat transfer problems at low temperature [66] or for low pressure gas [67] where heat conduction occurs mainly through phonon. Although for this particular cases the so-called Telegraphist equation is required, it usually brings unnecessary complexity for most solid and liquid, however in the case of of fluid of low thermal conductivity embedded in a high conductivity solid network, thermal transport cannot be described using a one-temperature model and LTE assumption.

If one considers a representative cell of such network called Control Volume (CV), and apply heat of one of its sides, it appears that the heat flux cannot be established instantaneously as described by the Fourier law. Instead, the heat first follows the strut of the network and then flows into the PCM. The mathematical formulation of this phenomenon was described by Tzou et al. [68] who considered such composite cell at a position r in space on which is applied a heat flux at a time t + τ q . Rather than following a straight line as they would do in a solid for example, the heat carriers first follows a longer conducting path leading to a delayed establishment of the heat gradient at a time t + τ T . In such case, the Fourier law becomes:

-k∇T(r, t + τ T ) = q(r, t + τ q ). (1.39)

Using a first order Taylor expansion, it leads to:

q(r, t ) + τ q ∂q ∂t r,t = -k∇T(r, t ) -kτ T ∂ ∂t ∇T(r, t ), (1.40)
in which both τ q and τ T are regarded as intrinsic properties of the composite. Substituting equation 1.40 into the heat equation gives the Dual Phase Lag Model [68]:

∂T ∂t + τ q ∂ 2 T ∂t 2 r,t = -α∇ 2 T(r, t ) -ατ T ∂ ∂t ∇ 2 T(r, t ), (1.41)
which introduced the dependence on the variation of the thermal gradient in time of the temperature field.

To give physical meaning to τ q and τ T , equation 1.41 can be obtained from the socalled Two-temperature model [19-21; 69; 70]

γ s ∂T s ∂t = ks ∇ 2 T s -h(T s -T f ), γ f ∂T f ∂t = k f ∇ 2 T f + h(T s -T f ), (1.42) in which γ s = (1 -ε)ρ s C p,s , γ f = ερ f C p, f , ks = (1 -ε)k s , k f = εk f and h is an integral heat
transfer coefficient which value is assumed to be independent of time but dependent of the thermal conductivities of both phases, on the heat transfer surface area, on the fluid volume ration (usually called porosity) and on the tortuousity of the solid network [START_REF] Vadasz | Lack of oscillations in Dual-Phase-Lagging heat conduction for a porous slab subject to imposed heat flux and temperature[END_REF].

This last coefficient has no physical sens but is rather a way to account for the change of scale: at the microscopic scale conduction is the actual way heat is transferred but looking at a macroscopic scale, two temperatures have to be considered and h is a mathematical way to link them while keeping as much information on the transfer as possible. If the LTE assumption applies then T f = T s and equation 1.42 is similar to equation 1.34, in the other case LTNE applies and T i (∀i ∈ [s, f ]) can be written:

T i = γ j h ∂T j ∂t - k j h ∇ 2 T j + T j , (1.43)
which can be substituted in the second term of equation 1.42 leading to:

α e ∇ 2 T i + α e τ T ∂ ∂t ∇ 2 T i -α e β e ∇ 4 T i = ∂T i ∂t + τ q ∂ 2 T i ∂t 2 , (1.44) with:
•

α e = k f +k s γ s +γ f , • τ q = γ f γ s h(γ s +γ f ) , • τ T = γ s k f +γ f k s h(k f +k s ) ,
•

β e = k f k s h(k s +k f ) .
Equation 1.44 is similar to equation 1.41 as long as β e can be neglected which is usually the case [START_REF] Vadasz | A paradox of heat conduction in porous media subject to lack of local thermal equilibrium[END_REF]. It appears that τ T and τ q depend on the two materials properties as well as on the heat exchange occurring at their interface and both vanish considering a perfect interface at which convective exchange occurs at an infinite rate. On the other hand, if the same coefficient vanishes (for example if the conductive matrix is coated with a highly insulating material), equation 1.42 cannot be coupled and the temperature in each phase is independent of the one in the other [START_REF]On the paradox of heat conduction in porous media subject to lack of local thermal equilibrium[END_REF].

CHAPTER 1. PHASE CHANGE MATERIALS

To consider LTNE, although it avoids any information loss in the homogenization process, adds complexity to the simulation and should only be used when necessary. Quintard et al. [20] determined that the validity of LTE was dominated by the quantity :

R = 〈T f 〉 f -〈T s 〉 s ∆〈T〉 , (1.45)
in which 〈T i 〉 i , i = s, f represents the intrinsic average temperature defined as:

〈T i 〉 i = 1 V V i T i d V, (1.46)
with V the total volume of the composite and V i the volume of a given phase and:

〈T〉 = ε f 〈T f 〉 f + ε s 〈T s 〉 s , (1.47)
with ε i the respective volume ratio of each phase. 〈T〉 represents the spatial average temperature of the composite. The term defined in equation 1.45 is the ratio of the temperature difference between the two phases at a given point in space and time to the variation of the average temperature on a characteristic length L T . If this quantity is sufficiently small (i.e. R << 1) then LTE is valid and expressing R is key. Quintard et al. [20] were able to express it as the ratio of two characteristic length: the first one l s f is characteristic of the heat exchange between the two phases while L is characteristic to the system size.

Equation 1.45 can be written as :

R = l s f L = 1 L × ε s ε f (ε f k s + ε s k f ) a v , (1.48)
with a convection-like heat exchange coefficient and a v the specific surface area (i.e.

the contact surface area between the fluid and the solid per volume unit). Both are linked to the previously described h as h = a v . To express h, Quintard et al. solved a so-called Closure problem on a particular case in which the solid structure is constituted of an array of cylinder. In this case: promote LTNE which could be explained if one looks at a representative volume of the cylinder array. Here, a cylinder of radius l s is enclosed in a cube of side length L. In this situation, the characteristic length of the fluid is expressed as l f = L/2-l s and the specific area can be defined as:

h = 8πε 2 f k f l 2 f ε s ε f k f ε s ks + 1 -3 ε f -4 ln ε s . ( 1
a v = 2πl s L 2 = π L 1 -2 l f L , (1.50)
which decreases with an increased l f /L ratio, limiting the contact surface area between the two phases hence the heat flow, delaying even more the establishment of the heat flux in the structure. Finally, the temperature difference seems to be sensible to the porosity as R tends to increase with it. This could be explained with the same example as before:

increasing porosity means decreasing the radius of the cylinder limiting again the contact surface area between both phases.

Hu et al. [58] simulated the melting of a PCM infused in a metallic foam under different conditions of power density using one and two-temperature models and compared the results to DNS.

Fusion-induced Convective Effects

Until now, only conduction has been seen as a heat transfer mode in the PCM modeling, but, as its temperature rises the PCM tends to expand, especially at the moment of the phase change (about 10 % in volume). As the PCM expands, buoyancy effects tend to appear, leading to convection in the melting material. To account for such phenomena, that will affect temperature distribution and melting rate, Voller [18] proposed a method considering three sets of equations dealing with the conservation of mass, momentum and energy. Before detailing the different sets of equation, one must define the superficial velocity (i.e. the average velocity) u:

u = f l u l , (1.51)
in which u l is the actual fluid velocity and f l is the liquid fraction defined in equation 1.27.

Then the momentum can first be defined as:

ρ ∂u ∂t + ρ(u • ∇)u = -∇P + µ∇ 2 u + S b + S m , (1.52)
with µ the fluid viscosity and P the pressure. S b is a source term used to depict the influence of buoyancy on the fluid motion. To do so, the Boussinesq treatment [START_REF] Nield | Convection in porous media[END_REF] is assumed to be valid [18]. This is true provided that density changes between the liquid and the the solid remain small in comparison with the actual material density throughout the flow region and that temperature variations are insufficient to cause the various properties of the medium to vary significantly. This supposes that the density of the melting material is constant everywhere except in this particular term in which the thermal expansion β is introduced leading to:

S b = ρg β(T -T m ), (1.53)
with g the gravitational force equivalent. On the other hand, the term S m account for the easier fluid motion as the fusion occurs. When the material is completely solid, it is easy to understand that no fluid motion can happen but as the temperature increases a mushy region (i.e. 0 < f l < 1) develops in which liquid and solid coexist. In this region, fluid motion is enabled and eased as melting occurs. This is usually modeled using the Darcy law that defines the superficial velocity as a function of the permeability K, depending on the liquid fraction f l :

u = - K µ ∇P, (1.54)
using this law, for low value of the liquid fraction, the permeability tends to zero so does the superficial velocity. In the numerical model, this behavior is accounted for by stating that:

S m = Au, (1.55)

CHAPTER 1. PHASE CHANGE MATERIALS

where A has a high value in the solid phase but tends to zero as the liquid fraction tends to one. Substituting this expression in equation 1.52, this source term takes a near zero value in the liquid region, and the momentum equation is in term of the actual fluid velocity. In the mushy region, however, A starts to increase such as the source dominates any other term in the equation and the velocity tends to zero as the liquid fraction approaches zero, which approximates the behavior of the Darcy law. The actual value of A was determined using the Carman-Koseny equation leading to :

A = -C (1 -f l ) 2 f 3 l + ω , ( 1.56)
in which C and ω are arbitrary big and small values, the former mainly depending of the morphology of the melting front (depends on the boundary conditions mostly) while the latter is introduced to avoid division by zero. Choosing those values, one should make sure that they allow for significant flow in the mushy region at high value of liquid fraction whereas A is large enough to suppress fluid velocities for near solid region.

The second equation of the set is called the Mass Conversion equation which simply states that no matter should escape the domain Ω on which the simulation is performed.

It is defined as follows:

∇ • u = 0. (1.57)
Finally, the last equation of the set is the Energy conservation equation that has been defined in the previous part considering LTE (see equation 1.34) or LTNE (see equation 1.42).

When fusion-induced convection is considered, these equations have to be modified to include the motion of the fluid. For example for LTE assumption, the Energy conservation equation becomes:

ρ eff C p,eff ∂T ∂t + ρ f C p, f u • ∇T = ∇ • (k eff ∇T). (1.58)
Considering fusion-induced convection modifies both temperature distribution and obviously melting rate and distribution. This has been shown in a number of cases [14; 59; 73; 75]. A good example of such modification is given by Yang et al. [START_REF] Yang | Comparison of direct numerical simulation with volume-averaged method on composite phase change materials for thermal energy storage[END_REF] who studied the melting of a PCM embedded in a copper foam. For that purpose, a representative model of such composite material is represented as a combination of several tetrakaidecahedron cell [START_REF] Thomson | On the Division of Space with Minimum Partitional Area[END_REF] (see figure 1.8) which are considered to be a good representation of an idealized foam cells. On this system, insulating boundary condition is applied on all faces except one, on which the temperature is set above the melting temperature of the PCM.

Results of the simulation are shown on figure 1.10 and 1.11, respectively in terms of melting fraction and temperature at different Fourier Numbers (F o = αt /L 2 ) for both DNS and volume averaged method (or homogenization method) considering LTNE assumption. Qualitatively, it can be said that in both cases the homogenization method is able to Figure 1.8: Computational domain and representative mesh for DNS [START_REF] Yang | Comparison of direct numerical simulation with volume-averaged method on composite phase change materials for thermal energy storage[END_REF]. yield similar prediction of the transient solid-liquid interface evolution although it lacks the pore-scale features the DNS can capture. For example, part of the PCM located close to the foam strut melts faster than the ones located in the middle of the cell: this cannot be captured by the volume averaging method that predicts a continuous melting front. Comparing now the convection to the conduction dominated simulation, it is observed that while the conduction dominated simulation predicts a melting front parallel to the heated wall, the shape predicted by the convection dominated simulation is different. Quantitatively, figure 1.9 depicts the total melted fraction for the four different cases mentioned before: it shows that both volume averaged method (VAM on the figure) and DNS convection dominated simulations produce higher melting rates (up to 40%) which is usually associated with a more homogeneous temperature distribution. This fact is reported in a number of studies [14; 59; 75; 77; 78] and is usually attributed to a better heat transfer due to the combined effect of conduction and convection.

A number of parameters can influence the establishment of convection in a PCM based system. It has been reported [62; 75] that foams, used as filler enhancing PCM thermal conductivity can significantly reduce natural convection due to viscosity and capillarity effects. This is especially true for foams composed of relatively small pores in which the space available for fluid motion is greatly reduced. Another parameter influencing the convection of the melted PCM is the value [59] and the direction [75; 79] of the gravity term relatively to the orientation of the PCM device. The influence of the value of the gravity term is easily understood as this is the driving force of convection. The orientation of the device compared to the gravity term affects the direction in which the fluid will preferably flow, if the direction is the same as the one promoted by heat conduction, the effect will combine but the further it is from the conduction direction the less help it brings to the melting process. This is shown on figure 1.12: a PCM enclosed in a cuboid box is heated from below. The experiment is repeated for different inclinations of the box and the shape of the melting front as well as the value of the melted fraction (see figure 1.13) are compared. The experiment shows that for increasing inclination angle the melting rate decreases.

Phase Change Materials Enhancement 1.4.1 PCM Fillers Types and Compared Enhancement Capacity

In the previous paragraphs, it has been repeatedly stated that for PCM to be used as LHT-ESS, they require to be enhanced in terms of thermal conductivity. To compare the different solutions and select the best one for a given application, two parameters need to be defined. The first one is called the effective thermal conductivity (ETC) which is the parameter used in the previous paragraph to represent the thermal conductivity of the composite formed by the filler and the PCM. It depends mostly on the thermal conductivity of the PCM and the base material forming the filler, the volume ratio of PCM inside the composite (called porosity ε) and the topology of the filler. Different models can be found in the literature depending on the filler and that question is discussed further for foams in the next paragraph. For now, the ETC tensor (it can be anisotropic) can be expressed as:

K eff =     f 1 (k s , k f , ε, G, ...), 0 0 0 f 2 (k s , k f , ε, G, ...) 0 0 0 f 3 (k s , k f , ε, G, ...)     , (1.59)
with f i ∈ [1,3] independent functions and G a parameter depending on the topology of the filler. As pointed out by Ji et al. [5] the ETC itself does not provide a full description of the filler performance. Although high thermal conductivity filler can enhance the PCM CHAPTER 1. PHASE CHANGE MATERIALS thermal conductivity, they reduce the amount of PCM composing the composite hence reducing its overall ability to store heat. A better, more complete parameter is referred to as the specific ETC (ϕ) and is expressed as:

ϕ = tr(K eff ) 3(1 -ε)k f , (1.60)
and can be seen as enhancing power of the filler normalized by its volume fraction.

The different solutions mentioned here can be divided into two groups: continuous and dispersed fillers. While dispersed fillers are composed of independent particles being randomly spread inside the PCM, continuous fillers usually form an actual network. Heat transfer wise, continuous fillers are expected to outperform dispersed ones as heat can flow from one side of the structure to the other without crossing any low conductivity zone whereas to go from one conductive particle or fiber of a dispersed filler, heat has to flow through PCM which causes additional thermal resistance and reduces the ETC for a given porosity. This is visible on figure 1.14 where Ji et al. [5] compared the specific ETC of continuous and dispersed fillers. It shows that for dispersed filler to reach the performance of continuous one, highly conductive materials like graphite or graphene (k

s 1000 W•K -1 •m -1 )
need to be used. On the other hand, it shows that, despite their relatively low thermal TPMS [START_REF] Schaedler | Architected Cellular Materials[END_REF] stituents but also by the spatial configuration [START_REF] Schaedler | Architected Cellular Materials[END_REF]. In contrast with foams in which spatial configuration is random, architected cellular materials can be tuned to optimize thermal or mechanical properties for example. This can be used for example to direct the heat flux in a given direction [9; 90], promote conduction over convection (or the exact opposite) [START_REF] Catchpole-Smith | Thermal conductivity of TPMS lattice structures manufactured via laser powder bed fusion[END_REF] or to increase the ETC for a given porosity compared to foams [9; 28; 86].

This last point is key in the development of efficient PCM fillers. In that sens, Ahmed outperformed by TPMS as it was the structure leading to the shortest time required for all the PCM to be melted. Ranking the different TPMS in terms of melting time showed that the IWP lead to the smallest time, followed by the Primitive and the Gyroid which was not expected due to the difference of ETC between the structures. The decorrelation between the values of ETC and melting time is related to the higher integral heat transfer coefficient of the IWP structure compared to the Primitive.

We showed here that a number of parameters can affect the overall performance of a given filler: nature (dispersed or continuous), base material (metal or carbon-based), solid volume fraction, topology or cell size. Although those proved to be useful in compar-

CHAPTER 1. PHASE CHANGE MATERIALS

ing the different fillers, they all required time consuming measurement and simulation.

To avoid it, analytical models are needed and would allow for an easier and quicker estimation of a given filler quality as well as for quicker simulation using the Homogenization theory (see part 1.3.2).

Effective Thermal Conductivity Models for Foams

To save time in both the assessment of a filler performance and in the simulation of the composite created, it is of interest to develop models able to predict the ETC of a given filler depending on its characteristic. Historically, and before the development of architected enhancers, foams have been widely used as filler due to their high specific area and convective heat transfer coefficient. Hence, the modeling of such filler has been of interest in the past decades and a number of models can be found in the literature.

Asymptotic models

The first type of models found in the literature are the asymptotic models [START_REF] Ranut | On the effective thermal conductivity of aluminum metal foams: Review and improvement of the available empirical and analytical models[END_REF] which provide two limiting conditions for the estimation of the ETC of a composite made of PCM and a solid conductive structure. The first one considers such material made out of layers of both materials oriented perpendicularly to the input heat flux (see figure 1.17). This model, called the Series model is expressed as:

k ⊥ = ε k f + 1 -ε k s -1
.

(1.61)

On the other hand, using the Parallel model, the layers are supposed to be oriented in the direction of the heat flux leading to: Although relatively simple, those cannot account for the conductivity of foam as its strut are usually oriented in random directions of space. As a consequence equations 1.61 and 1.62 are usually considered to be respectively a lower and an upper bond of foam ETC.

k ∥ = k f ε + k s (1 -ε). (1.62) (a) (b) (c) (d)
One way to use those model anyway is to express the ETC as a combination of both, adding an empirical variable F accounting for the fraction of solid oriented in one direction or another. In this situation k eff is written:

k eff = k F ∥ k (1-F) ⊥ , 0 ≤ F ≤ 1. (1.63)
Finding the value of F is done solving the following equation:

F = ln ε k eff k f + (1 -ε) k eff k s ln 1 + ε(1 -ε) k s k f + k f k s -2 , ( 1.64)
which was solved by Singh and Kasana [START_REF] Singh | Computational aspects of effective thermal conductivity of highly porous metal foams[END_REF] and improved by Kumar and Topin [START_REF] Kumar | Determination of effective thermal conductivity from geometrical properties: Application to open cell foams[END_REF] using experimental data, leading to :

F = 0.3031 + 0.0623 ln ε k s k f . (1.65)
Using this model and a set of 2000 experimental values, Kumar and Topin [START_REF] Kumar | Determination of effective thermal conductivity from geometrical properties: Application to open cell foams[END_REF] could predict the ETC of high porosity foam (ε > 0.9) in the error range ±6%. Although those models are efficient predicting the ETC of foam, they require a great number of experimental data and could be biased if the set is not well chosen. For that reason, we discuss the possibility of using analytical models.

Analytical models

Following the growing development of foam as PCM filler, an increasing number of analytical models could be found in the literature. All of this models follow the same process:

the foam is considered to be the repetition of a unique ideal cell on which a given heat flux q is applied on one of its sides. Then, the cell is sliced along the heat flux axis in layers on which the ETC can be determined using the Parallel model (eq. 1.62). Finally, the Following this process, Calmidi and Mahajan first determined the solid (V s ) and fluid (V f ) volume fraction in the first layer, leading to:

V s = t (L + b)w, V f = t (3L -(L + b)) w.
(1.66)

with w the hypothetical width in the direction perpendicular to the plan. Substituting CHAPTER 1. PHASE CHANGE MATERIALS equation 1.66 into equation 1.62, the ETC of the first layer is expressed as:

k I = k f + k s -k f 3 1 + b L . (1.67)
Similarly the ETC of layers 2 and 3 are expressed as follows:

k II = k f + 2 3 (k s -k f ) b L , k III = k f + 4 3 3 (k s -k f ) t L .
(1.68)

The individual ETC of each layer are combined using the Series model, hence:

L 1 + L 2 + L 3 k eff = L 1 k I + L 2 k II + L 3 k III . (1.69)
which is equivalent to:

k eff = 3 2   e α k f + k s -k f 3 (1 + e) + (1 -1/α)e k f + 2 3 (k s -k f )e + 3 2 -e k f + 4/α 3 3 (k s -k f )e   -1
.

(1.70)

In equation 1.70, α = b/t is strut to node cross sectional area ratio that is constant with porosity. On the other hand, e = b/L is the node to strut length ratio which is key in determining the ETC of the foam but is also difficult, if not impossible, to measure. To determine the value of both e and α, Calmidi and Mahajan first expressed e as a function of α and the porosity ε that can be written:

ε = 1 - 2 3 t (L + b) + 2b(b -t ) + ( 3 2 L -t ) 4t 3 3L 2 , ( 1.71)
hence:

e = b L = -1 α + 1 α 2 + 2 3 (1 -ε) 2 -1 α 1 + 4 3 2 3 2 -1 α 1 + 4 3 (1.72)
Comparing the previously described equations and experimental data, Calmidi and Mahajan determined that a ratio α = 1/0.09 resulted in excellent agreement even though it was noted that this was only a representative value that usually varies depending on the node one observes.

Going further into the development of precise, reliable ETC model for foams Bhattacharya et al. [43] improved Calmidi and Mahajan model by considering a similar representative network but, this time, considering circular nodes instead of square ones. This modification requires to modify the way the unit cell is sliced as shown on figure 1.20.

Similarly, the thermal conductivity of each layer can be calculated leading to:

k I = k f + k s -k f 3Lt   Lt + R 2 2   sin -1 t R + t R 1 - t 2 R 2     , (1.73) CHAPTER 1. PHASE CHANGE MATERIALS Layer III Layer I 2t Layer II Layer IV Layer V t R L cos60 3L 60� Figure 1
.20: Bhattacharya representative medium unit cell [43].

k II = k f + R 2 (k s -k f ) 3L 2 3(R 2 -t 2 ) -3t × sin -1 1 2R 3(R 2 -t 2 ) -t -sin -1 t R + 1 2 sin 2 sin -1 1 2R 3(R 2 -t 2 ) -t - 1 2 sin 2 sin -1 t R , (1.74
)

k III = k f + k s -k f 3Lt × Rt + R 2 2 sin -1 y 2 R -sin -1 y 1 R + 1 2 sin -1 2 sin -1 y 2 R - 1 2 sin -1 2 sin -1 y 2 R + 3 2 y 2 2 -y 2 1 2 - 3 2 y 1 t , (1.75
)

k IV = k f + k s -k f 3L(R -y 2 ) R 2 π 2 -sin -1 y 2 R - R 2 2 sin 2 sin -1 y 2 R + 4t 3 (R -y 2 ) , (1.76
)

k V = k f + (k s -k f ) 4t 3 3L . (1.77)
with R the strut radius,

y 1 = 1 2 ( 3(R 2 -t 2 )-t ) and y 2 = 1 2 ( 3(R 2 -t 2 )+t ).
Once again, the determination of α = R/t is key to precisely predict the behavior of the foam. Rather than fitting the model to experimental values, Bhattacharya et al. chose to rely on physical measurement from pictures taken under microscope and determined that 1/α was in the range 0.17 to 0.21 with an average value around 0.19. Comparison with experimental results showed once again good agreement.

Although the last two models described here showed good agreements with experimental values, they are based on a two-dimensional representation of the foam. To get closer to the actual structure of the foam which would improve model accuracy, a three dimensional representation of the foam is necessary. Based on work from Lord Kelvin [START_REF] Thomson | On the Division of Space with Minimum Partitional Area[END_REF], it was determined that the Tetrakaidecahedron (see figure 1 and eight hexagons and is the idealized shape that will most likely be attained by the foam from the nature of the foaming process. Experimental work based on the observation of soap foam by Matzke [START_REF] Matzke | The Three-Dimensional Shape of Bubbles in Foam-An Analysis of the Role of Surface Forces in Three-Dimensional Cell Shape Determination[END_REF] tended to confirm Kelvin's hypothesis. Kelvin's work was more recently corrected by Weaire and Phelan [START_REF] Weaire | A counter-example to kelvin's conjecture on minimal surfaces[END_REF] (see figure 1.21b) and then Gabbrielli [START_REF] Gabbrielli | A new counter-example to Kelvin's conjecture on minimal surfaces[END_REF] (see figure 1.21c) but the relative simplicity of the Kelvin cell combined with the relatively small differences in terms of specific area and volume of the other cells made it the best compromise for the calculation of the ETC.

Based on Kelvin's assumptions, Boomsma and Poulikakos [START_REF] Poulikakos | On the effective thermal conductivity of a threedimensionally structured fluid-saturated metal foam[END_REF]

V A,s = e 2 + 1 2 d π(1 -e) d L 3 , V B,s = 1 2 e -d e 2 L 3 , V C,s = 1 -e 2 πd 2 L 3 , V D,s = 1 4 e 3 L 3 , (1.78)
accounting for the corrections by Dai et al. [44].

Using the previous equation, the porosity can be defined as:

ε = 1 - 2 2 d e 2 + πd 2 2 (1 -e) + e 2 -d e 2 + πd 2 (1 -e 2) + e 3 4 .
(1.79)

As shown in equation 1.59 the final goal is to express the ETC as a function of the thermal properties of the PCM and the conductive material and of the porosity, hence d and e need to be eliminated in the final expression. To do so, Boomsma and Poulikakos expressed d as follows:

d =   2(2 -3 2 4 e 3 -2ε) π(3 -4e 2 -e)   1 2
.

(1.80)

Then, equation 1.69 is re-written in term of thermal resistance R i , i ∈[A,B,C,D] of each layer, leading to:

k eff = 2 2 R i , (1.81)
with

R A = 4d (2e 2 + πd (1 -e))k s + (4 -2e 2 -πd (1 -e))k f , R B = (e -2d ) 2 (e -2d )e 2 k s + (2e -4d -(e -2d )e 2 )k f , R C = 2( 2 -2e) πd 2 2k s + 2(2 -πd 2 2)k f , R D = 2e e 2 k s + (4 -e 2 )k f , (1.82)
accounting for the corrections by Dai et al. [44].

Before the ETC can actually be calculated, the value of e needs to be estimated. To do so, Dai et al., who corrected and extended Boomsma and Poulikakos model, once again fitted their model against experimental data and found that e = 0.198 was the best fit.

The last model of interest found in the literature was developed by Yang et al. [7].

It uses a similar unit cell and slicing method although the calculation method is slightly

Q t = k eff A 0 ∆T H c , (1.83)
with A 0 the area of the unit cell top surface and ∆T the temperature difference between the top and the bottom of the cell. Then, Yang considered a parallel heat conduction between the solid ligaments (Q s ) and the fluid (Q f ) which can be mathematically repre-

sented by Q t = Q s + Q f .
Due to the hypothetical high porosity of the foam Q f can simply be written:

Q f = k f A f ∆T H c , (1.84)
with A f the fluid cross sectional area along the heat flux axis, considered constant along the cell. On the other hand, Q s expression is more complicated due to the intricate network formed by the foam, hence the use of an integral form is necessary:

Q s = -k s A s (s) d T d s , (1.85)
where s is the axis along the foam strut depicted on figure 1.23 and A s (s) the solid cross sectional area at point s. Using the separation of variables method the previous equation becomes:

Q s A s (s) d s = -k s d T, (1.86)
which is integrated along the cell, leading to:

Q s = -k s ∆T/ H c 0 1 A s (s) d s. (1.87)
Combining the respective expression of Q t , Q f and Q s , it leads to:

k eff = H c A 0 H c 0 1 A s (s) d s -1 k s + εk f , (1.88)
if the approximation is made that A f /A 0 = ε. To go further, it is necessary to estimate the integral in the previous equation. For that purpose two parameters are first defined:

e = t /L, e ≥ 0, α = A t /A L , α ≥ 1.
(1.89)

CHAPTER 1. PHASE CHANGE MATERIALS

as well as the porosity: ). Using this, the integral can be calculated as:

ε = 1 - A L L 2 × 6(1 -e) + 3αe 4
H c 0 1 A s (s) d s = 3 i =1 I i = 3 i =1 z i z i -1 1 A s (s) d s (1.91)
with z 0 = 0 and z 3 = 2L/2. Then:

I 1 = t /2 0 1 A s (s) d s = t A t = eL αA L , I 2 = 2L/2-t /2 t /2 1 A s (s) d s = L -t A L = L(1 -e) A L , I 3 = 2L/2 2L/2-t /2 1 A s (s) d s = t 2A t = eL 2αA L , (1.92)
and finally :

H c 0 1 A s (s) d s = L A L × 1 -e + 3e 2α . ( 1
k eff = 1 -ε 1 -e + 3e 2α 3(1 -e) + 3 2 eα k s + εk f . (1.94)
Once again, the key point in assessing the value of the composite ETC lies in the correct determination of both α and e. Yang et al. determined, comparing them with experimental values, that α = 3/2 and e = 0.3 lead to satisfying results as the Root-Mean-Square (RMS) deviations from the experimental values lied between 9.8% and 11.1% for foam of porosity ε ≥ 0.9.

All previously described models consider constant values for both node to strut cross sectional area ratio (α) and for the node to strut length ratio (e). Although this was done to fit experimental data, this is physically impossible. This issue was first discussed by Dai et al. [44], and then Yang et al. [7] who both pointed out that in any foam the edge of the node has to be be greater than the strut diameter leading to t /a ≥ • First: Yang et al. [7] with a RMS deviation of 11.1%,

e = a + bε + cε 2 + d ε 3 , (1.95)
• Second: Boomsma and Poulikakos [START_REF] Poulikakos | On the effective thermal conductivity of a threedimensionally structured fluid-saturated metal foam[END_REF] (corrected by Dai et al. [44]) with a RMS deviation of 12.05%,

• Third: The model developped by Sing and Kasana [START_REF] Singh | Computational aspects of effective thermal conductivity of highly porous metal foams[END_REF] based on the combination of the Series and Parallel models with a RMS deviation of 13.96%,

• Fourth: Calmidi and Mahajan [START_REF] Calmidi | The Effective Thermal Conductivity of High Porosity Fibrous Metal Foams[END_REF] with a RMS deviation of 20.3%.

From there, it seems that models based on a 3D unit-cell outperform the ones based on a 2D network justifying the need for better ways to represent the structure of the foam to get the best estimation of its ETC. .

Numerical Models

Following the last paragraph's conclusion, it is clear that a better representation of the foam would lead to a better estimation of the ETC as well as a better understanding of the foam's features influencing it. Previously, foams were represented by a unique and ideal unit cell : the Tetrakaidecahedron. Although, this representation has been proven to approach the foam average unit-cell and to predict quite precisely its ETC, it prevents one from accounting for the local features of the foam and the actual shape on its struts and nodes.

In the previous part, the foam has been represented by an ideal unit-cell with cylindrical strut and cubic nodes. In addition, the cell size is unique and the distribution of cell size inside the foam is not accounted for. However, an actual foam is different from its ideal representation as it was pointed out by Jang et al. ing on the minimum to maximum cross sectional area ratio on a given strut as well as the associated porosity that remains relatively constant. Finally, the influence of the average cell size and its distribution is discussed. It was shown by most studies that the average cell size does not affect the ETC of foam no matter the material, porosity or strut profile.

However, Baillis et al. [45] studied the influence of the cell size distribution. As mentioned earlier, the cell size across the foam follows a log-normal law, with a mean value Ps which corresponds to the average cell size and a given standard deviation σ. Using those two values, Baillis et al. defines the normalized standard deviation σ n as:

σ n = σ Ps . ( 1
k eff = (0.9 -0.33 × σ n ) × (1 -ε) 2 + 1 3 (1 -ε) k s . (1.97)
This graph clearly shows the dependency of the foam ETC on the normalized standard deviation with a decreased value for more poly-dispersed structure, with the highest values obtained for ideal Tetrakaidecahedron structure.

While analytical models were already able to predict the dependency of the foam ETC on the porosity, numerical models are able to account for the statistical features of the foam like the strut profile or cell size distribution. Those model showed that higher ETC was obtained for cylindrical strut with minimum matter concentration at the nodes and for less poly-dispersed networks, with the best results obtained for periodic ones.

System Level Enhancement : Cascaded Storage

Previously, LHTESS systems were described as the combination of a PCM for latent heat energy storage and a filler for thermal conductivity enhancement. The filler in which the PCM is embedded can be of different nature : nano-particles, foam, fin, pin, with varying size and conductive material volume fraction, each solution providing different advantages and drawbacks. One way to overcome the drawback brought by one solution may be to combine several ETC enhancement solutions to produce a more complex and better filler. Such solutions are found in the literature in the form of combined fin and foam enhancers [12; 14], foam displaying porosity gradient in one or two directions or material variations [13; 106-108] or a combination of foam and nano-particles [START_REF] Mahdi | Solidification enhancement with multiple PCMs , cascaded metal foam and nanoparticles in the shell-and-tube energy storage system[END_REF].

All these solutions have been proven to enhance both charging/discharging (melting/solidifying) time as well as promoting a more homogeneous temperature in the system. This was clearly shown by Zheng et al. [START_REF] Zheng | Effect of metal foam with twodimensional porosity gradient on melting behavior in a rectangular cavity[END_REF] who studied the influence of foam porosity gradient in a two-dimensional square enclosure heated up on the left side. tion dominated with the gravity oriented vertically. In this case, the horizontal gradient only brought small melting time reduction as the melting front speed increases close to the heat source while it is reduced on the opposite wall. On the other hand, the vertical gradient allows an easier liquid convection at the top of the enclosure where it usually happens while it brings higher ETC where its effects is less felt. The two dimensional gradient, finally, combined both effects which produced the lowest melting time.

On the other hand, Yang et al. The relatively thick and conductive pins are rapidly cooled down creating, artificially increasing the contact surface area with the cold bottom surface, enabling a better heat conduction from the PCM where the heat is stored, to the heat sink.

All examples from the literature showed similar results: the performance of a system are not only linked to relative volume fraction of conductive material and PCM but also to the location of each one, relative to the heat source. It seems to the authors that conductive material should be placed close to the heat source, as way to rapidly spread the heat that can be stored later and further away from the source as latent heat.

CHAPTER 1. PHASE CHANGE MATERIALS Figure 1.30: Pin-foam enhanced PCM system [12].

Conclusion

In this first chapter, we draw the baseline of the study of PCM based system. First, we defined what was a PCM by defining its basic properties: adapted melting point, high latent heat of fusion, high specific heat, high thermal conductivity, high nucleation rate and crystal growth, reversible freezing/melting cycle, contain no flammable or poisonous compounds and finally, be available in large quantity and at low cost. Then, we described the wide variety of chemical compositions: from organic to metallic or other inorganic compounds. The change in chemical composition was obviously linked to a wide change in thermal properties, especially in terms of latent heat of fusion and melting temperature. Those two parameters can be wisely chosen to adapt to the wide variety of applications of PCM based systems: from home thermal comfort in which both heat flux applied to the system are mild, to hybrid or electrical battery or high power electronic thermal management with both higher heat flux and temperature. Through the study of these examples, three issues arose : PCM low thermal conductivity limiting their performance, fusion influence and modeling which is a highly non-linear phenomenon that can be accelerated through buoyancy, and finally the relation between the PCM and the filler aiming at increasing its thermal conductivity and how both can be wisely modeled and simulated to save calculation time.

The second part of the chapter was dedicated to the study to both fusion model and homogenization theory. We detailed here, the different methods aiming at modeling the fusion of a given material from the Stefan [17] method: a first simplified technique that has shown some success in the past to the fixed grid enthalpy method of Voller [40] that allows relatively simple and precise modeling of a melting material if its melting temperature and heat of fusion are known. Then, we discussed the Homogenization theory or the way the composite formed by a given PCM and a filler aiming at enhancing its thermal conductivity can be modeled as one and unique effective bulk material for time saving purpose. Another aspect of homogenization is the relation between the PCM and the filler which topology (mainly specific area) can influence the heat spreading dynamic and increase the melting rate of the PCM if not adapted. Finally, we depicted the influence of buoyancy in the melting process as the denser solid PCM melts, the less-dense liquid PCM shows buoyancy induced convection that usually tends to speed up the melting process as it helps homogenizing the temperature in the system.

The third and last part of the chapter intends to better understand what should be the characteristic of a suitable filler. For that purpose, we first compared the discrete to the continuous fillers and concluded that the discrete filler were outperformed by the continuous ones due to the lack of inter-connectivity of the conductive particles adding thermal CHAPTER 1. PHASE CHANGE MATERIALS resistance. Among the continuous fillers, it was shown that the topology of the representative cell of the medium played a undeniable role in the enhancement performance of the filler. It was proven, for example, that TPMS structures would outperform foam, that are the most used filler nowadays. To better understand this last fact, we described the different ways aiming at predicting the foam ETC, from the analytical one that give a quick and relatively precise prediction to the numerical ones that, while being more time consuming, allowed for refined interpretations on the necessity to use architected fillers. Finally, we showed that the choice of the optimum filler with the optimum material and porosity could be further enhanced by a smart positioning of conductive material (composing the filler) and PCM in order to better spread the heat from the heat source leading to a better homogenization of the temperature on the system level inducing a higher melting rate and better overall performances.

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

Introduction

As mentioned in the previous chapter, even-though foams possess a high specific area (area per unit volume), high porosity and can be made of highly conductive material making them great candidates as heat sinks or in our case Phase Change Material (PCM) enhancers. However, they suffer from their random structure and the way they are manufactured. Actually, the latter tends to create hyperbolic struts [1; 2] with reduced cross section in their middle being the cause of additional thermal resistance while the former generates a dispersion in terms of pore size distribution being the cause of a reduced Effective Thermal Conductivity (ETC) [3][4][5]. From those statements, it seems that architected or structured material could reach higher ETC as both strut profile, strut orientation, strut positioning and pore size distribution could be controlled.

Architected materials have seen a huge development these last few years due to their great properties. Actually, their properties depends both on the material they are made of, and on the spatial configuration of solid and void [6]. These materials have been widely used due to their outstanding combination of lightweight and mechanical properties especially in the aviation industry but also in the biomedical industry as scaffolds for tissue engineering as they have bone-like structure [7]. More recently, engineered or more generally cellular materials have been used because of their thermal properties. Depending on the intrinsic properties of the material and their topology, they can be used as thermal insulator, thermal enhancer in heat storage unit or heat exchanger. In regard of the previous paragraph, the use of such engineered materials could lead to a better thermal conductivity as they have a perfectly periodic structure avoiding any of the effects previously mentioned. A large variety of periodic cellular material exists: we can cite lattice structures [8][9][10] , shellular materials [11][12][13][14], hexagonal honeycomb, fins... Another reason of their development is linked to the recent progress in additive manufacturing which has allowed manufacturing of more complex (like Triply Periodic Minimal Surface (TPMS) [15], Lattice [9], or bio-inspired shapes [16]...) and smaller [17; 18] patterns made in a greater variety of materials [19][20][21].

Due to their relative simplicity and given the freedom they offer in terms of topology tuning, we here concentrate on lattice structures. In this chapter, an ETC model for lattice structures is developed in order to assess the possible improvement due to topology changes from foam, to reduce computing time during simulations using the Local Thermal Equilibrium (LTE) hypothesis described in the previous chapter and finally to optimize PCM based thermal management devices. To do so, we will first present the mathematical definition of a model of ETC for some lattices, inspired by the work of Yang et al. [22] on foams. In this first part, porous network are modeled as isolated single unit cells which allow to compare them but conceal any side effect that may occur. In a sec-

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

ond part, a model of thermal constriction, being the cause of these side effects is implemented into the unit cell model to account for the modification of the ETC when more realistic structures are dealt with. Following this, a series of simulations aiming at theoretically validating the previous model is detailed. Finally, a new method to process flash laser method measurement results is developed to measure the thermal conductivity of a series of samples, specifically designed to validate the previously developed models.

Effective Thermal Conductivity Model

Lattice Structures Effective Thermal Conductivity

As mentioned in the previous paragraph, some foam features reflecting its randomness lead to a reduction of the ETC. It seems coherent to think that periodic structure with optimized topology would enable an ETC increase for a given porosity. To confirm this hypothesis four different lattices (BCC, FCC, BCC z and FCC z ) depicted in figure 2.1 are considered because of their relatively simple topology and the fact that they can easily be manufactured through additive manufacturing. It is to be noted that this model could be adapted to any other lattice structure. To calculate lattice structures ETC, the model developed by Yang et al. [22; 23] which is detailed in chapter 1, is adapted. That particular model was preferred to others due to its relative simplicity and flexibility regarding the structure to be studied. For clarity the first steps of the calculation are reminded in the following.

As it was done by Yang et al., we consider a single unit-cell of a larger structure to calculate the ETC. This assumption, that was considered correct for foam, seems coherent

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

anisotropic which involves the definition of two different ETC : k z in the vertical direction (or out-of-plane) and k x y in the horizontal plane (or in-plane direction). To determine their ETC we first consider a heat flux Q applied on the superior face (in red) of a given unit cell of side P s . Using Fourier law we can link the heat flux to the ETC k z and the temperature difference between the superior and the inferior faces ∆T using:

Q = k z A 0 ∆T P s . (2.1)
The amount of heat transferred to the CV is conducted along the solid (Q s ) and the fluid (Q f ) phase in parallel, leading to :

Q = Q s + Q f . (2.
2)

It should be noted that this assumption is verified only if k s k f . Along the CV length the cross sectional area of the fluid phase stays almost constant due to the relatively large fluid volume ratio, called porosity (ε). Hence we can write:

Q f = k f A f ∆T P s . ( 2.3)
On the other hand, due to its complex geometry the cross sectional area of the solid phase varies along the CV, here a global expression of the heat flux over the entire cell cannot be guaranteed and a differential form along the strut axis s is preferred:

Q s = -k s A s (s) d T d s . (2.4)
Applying the separation of variables method and integrating along the CV, yields:

Q s = k s ∆T/ P s 0 1 A s (s) d s. (2.5)
Finally we can substitute eq. (2.1), (2.3) and (2.5) in eq. (2.2) leading to:

k z = k s P s A 0 / P s 0 1 A s (s) d s + k f ε. (2.6)
For simplification purpose, eq. (2.6) is written as:

k z = k f ε + k s G s,z (1 -ε), (2.7)
with :

G s,z = 1 1 -ε P s A 0 / P s 0 1 A s (s) d s . (2.8)
In equation 2.8, we define G s,z (a similar G s,x y will be defined later) which is an image of the influence of the topology on the ETC. We actually see that equation 2.7 only differs from mixing law by this unique term that could be defined as a non-dimensional thermal resistance. To calculate it, the schemes depicted in figure 2.1 are used. On those ones, the four different lattices (i.e. BCC, FCC, BCC z and FCC z ) are represented as a combination of cylinders, representing their struts and cubes representing their nodes. This method is similar to the one Yang et al. [22] used for foams and has been proven reliable [22][23][24][25]. Following Yang's formalism, struts dimensions are described by a set of dimension (A l and L) respectively representing their cross sectional area and length while nodes are represented by a similar set (A t and t ) in which t denotes their height. It should be noted that due to their respective topology, all cells, except the BCC one, possess two sets of nodes; their respective dimensions are signaled by an index: 1 for the corner nodes or 2

for the centered or face centered nodes.

To simplify G z (and later G x y ) calculation, the different cells are sliced into N different parts figure 2.2); which provides a way to differentiate zone where nodes are present to the other where only struts can be found. Using the slicing, we are able to express:

P s 0 1 A s (s) d s = N i =1 S i 1 A s (s) d s.
(2.9)

Profiting from the cell symmetry we can write:

P s 0 1 A s (s) d s = 2 × t 1 /2 0 1 A s (s) d s + (P s -t 2 )/2 t 1 /2 1 A s (s) d s + P s /2 (P s -t 2 )/2 1 A s (s) d s . (2.10)
The way equation 2.9 is expressed, depends on the cell we consider. Actually depending on the cell type, struts and node arrangement will vary, modifying the expression of

P s 0 1 A s (s) d s.
In the following paragraphs we detail its calculation for all four cell in the outof-plane and in-plane direction as well as the expression of porosity as a function of their respective strut and nodes dimension.

BCC Cell: Porosity and Topology Parameters Calculation

Starting from equation 2.10, we express each of the three terms from the right side of the equation as functions of the struts and nodes dimensions. As BCC cell is an isotropic cell

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

that can be described using only one type of node, calculations are simplified and G s,z and G s,x y are identical. For this cell, we have:

t /2 0 1 A s (s) d s = t 2A t , (2.11)
(P s -t )/2 t /2 1 A s (s) d s = L -ω z t 4A l , (2.12)
P s /2 (P s -t )/2 1 A s (s) d s = t 2A t . (2.13)
In equation 2.12 an additional ω z term is added. This term, whose value will be determined in a following part, helps accounting for the complex interface between nodes and struts that make the actual struts and nodes lengths determination difficult. For anisotropic cells an additional ω x y will be added. To simplify the expression of the last three equations we define, as Yang et al. did, two sets of dimensionless parameters relating strut and node dimensions as follows:

α = A t /A l , α > 1, (2.14) e = t /L, e ≥ 0. (2.15)
As for ω z , the value of these parameters will be determined and discussed in the next part.

Substituting equations 2.11 to 2.15 into equation 2.10, we derive:

P s 0 1 A s (s) d s = L 2A l 1 -ω z e + 4e α . (2.16)
To be able to determine the complete expression of G s,z we now need to express the porosity as a function of the dimensionless parameter of the cell. Porosity is defined as the ratio of the void space (i.e. non-solid) in the cell to the cell total volume (V t ot ). To express it, we first determine the total volume of solid (V s ) in the cell:

V s = 8 × (L -ω ε t )A l + 2 × t A t , (2.17)
as well as the total volume V t ot :

V t ot = P 3 s = 8L 3 3 3 , ( 2.18)
with ω ε another scaling constant assigned to porosity calculation. This last equation can be written, using the dimensionless parameters, as:

V s = 8 × LA l (1 -ω ε e + eα 4 ). (2.19)
Then, we derive: for the BCC cell:

ε = V t ot -V s V t ot = 1 - A l L 2 3 × 3(1 -ω ε e) + 3
G s,z = 1 -ω z e + 4e α 3(1 -ω ε e) + 3 4 αe -1 . ( 2

FCC Cell: Porosity and Topology Parameters Calculation

We adopt a similar method for the FCC cell except that this time two parameters have to be calculated to account for the cell anisotropy. In addition, two types of nodes can be found in that cell requiring the use of two sets of dimensionless parameters (i.e. (α 1 , e 1 ) and (α 2 , e 2 )). Following the same method than previously we obtain, in the out-of-plane direction (see figure 2.3a):

t 1 /2 0 1 A s (s) d s = t 1 2A t ,1 , (2.22)
(P s -t 2 )/2

t 1 /2 1 A s (s) d s = L -ω z (t 1 + t 2 )/2 4A l , ( 2
P s 0 1 A s (s) d s = L 2A l 1 -ω z e 1 + e 2 2 + 2 × e 1 α 1 + e 2 2α 2 . (2.25)
In a similar way than previously, we need to determine the expression of porosity for the FCC cell. We first express the total volume of the cell as a function of L:

V t ot = P 3 s = 2 2L 3 . (2.26)
Then we derive the solid volume expression in such cell:

V s = 8 × L -ω ε t 1 + t 2 2 + t 1 A t ,1 + 2 × t 2 A t ,2 f (2.27)
from which we derive the expression of the porosity: We finally derive the expression of G s,z for a FCC cell:

ε = 1 - A l L 2 2 2 1 -ω ε e 1 + e 2 2 + e 1 α 1 + 2e 2 α 2 4 . ( 2
G s,z = 2 1 -ω ε e 1 + e 2 2 + α 1 e 1 + 2α 2 e 2 4 1 -ω z e 1 + e 2 2 + 2 e 1 α 1 + e 2 2α 2 -1
.

(2.29)

As mentioned previously, FCC cell is anisotropic which, contrary to the BCC cell, requires the calculation of an additional term G s,x y from which we can derive the in-plane ETC. For that purpose, we have to change the way we slice the cell as depicted on figure 2.3b. Using the same method than for the out-of-plane direction, we derive:

t 1 /2 0 1 A s (s) d s = t 1 2A t ,1 , (2.30)
(P s -t 2 )/2

t 1 /2 1 A s (s) d s = L -ω x y (t 1 + t 2 )/2 2A l , ( 2
P s 0 1 A s (s) d s = L A l 1 -ω x y e 1 + e 2 2 + e 1 α 1 + e 2 α 2 . (2.33)
The value of porosity remains the same than the one in equation 2.28, bringing:

G s,x y = 2 1 -ω ε e 1 + e 2 2 + α 1 e 1 + 2α 2 e 2 4 × 2 1 -ω x y e 1 + e 2 2 + e 1 α 1 + e 2 α 2 -1
.

(2.34)

BCC z Cell: Porosity and Topology Parameters Calculation

We follow the exact same process to derive G s,z and G s,x y for the BCC z cell. Once again we start with the evaluation of P s 0 1 A s (s) d s, in the out-of-plane direction (see slicing on figure 2.4a), based on equation 2.10:

t 1 /2 0 1 A s (s) d s = t 1 2A t ,1 , (2.35)

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

(P s -t 2 )/2

t 1 /2 1 A s (s) d s = L -ω z (t 1 + t 2 )/2 (4 + 3)A l , (2.36) P s /2 (P s -t 2 )/2 1 A s (s) d s = t 2 2(A t ,2 + A l ) , (2.37)
from which we derive:

P s 0 1 A s (s) d s = L A l 2 1 -ω z e 1 +e 2 2 4 + 3 + α 1 e 2 + (1 + α 2 )e 1 α 1 (1 + α 2 ) . (2.38)
Anew, we determine the solid volume in the cell:

V s = 8 × A l L -ω ε t 1 + t 2 2 + A l ( 2 3 L -t 1 ) + e 1 α 1 + e 2 α 2 , (2.39)
leading to:

ε = 1 -3 A l L 2 3 1 -ω ε e 1 +
G s,z = 2 × 3 1 -ω ε e 1 + e 2 2 + 3 4 + 3 8 (α 2 e 2 + e 1 (α 1 -1)) × 2 1 -ω z e 1 +e 2 2 4 + 3 + α 1 e 2 + (1 + α 2 )e 1 α 1 (1 + α 2 ) -1
.

(2.41)

To carry out G s,x y calculation for a BCC z cell, we need to change slicing strategy. Actually as shown in figure 2.4b, two additional slices were necessary to account for the fact that the supplementary struts in the out-of-plane direction may also participate to the heat conduction in the in-plane direction. To properly derive P s 0 1 A s (s) d s in that situation we need to define strut radius r = A l /π and pore size to strut radius ratio β = P s /r . Anew, using cell symmetry we can evaluate the latter integral by:

r 0 1 A s (s) d s = 2L/β 3 A t ,1 + 2 3 L -t 1 2r π , (2.42
) We finally carry out both G s,z and G s,x y calculation for the FCC z cell. The calculation is similar in all points with the one performed for BCC z cell. We detail here the same steps starting with P s 0 1 A s (s) d s evaluation in the out-of-plane direction:

t 1 /2 r 1 A s (s) d s = t 1 -4 3β L 2A t ,1 , (2.43)
(P s -t )/2 t /2 1 A s (s) d s = L -ω z t 4A l , ( 2
t 1 /2 0 1 A s (s) d s = t 1 2A t ,1 , (2.47)
(P s -t 2 )/2

t 1 /2 1 A s (s) d s = L -ω z (t 1 + t 2 )/2 (4 + 2)A l , (2.48)
P s /2 (P s -t 2 )/2 1 A s (s) d s = t 2 2(2A t ,2 + A l ) . (2.49)
Substituting in equation 2.10, we obtain:

P s 0 1 A s (s) d s = L A l 2 1 -ω z e 1 +e 2 2 4 + 2 + α 1 e 2 + (1 + 2α 2 )e 1 α 1 (1 + 2α 2 ) . ( 2

.50)

As it was previously done for other cells, we now define the solid volume inside a FCC z cell:

V s = 8 × A l L -ω ε t 1 + t 2 2 + A l 2L -t 1 + e 1 α 1 + 2 × e 2 α 2 , (2.51)
leading to the expression of the porosity in such cell:

ε = 1 -2 A l L 2 2 1 -ω ε e 1 + e 2 2 + 1 2 2 + 2α 2 e 2 + e 1 (α 1 -1) 4 .
(2.52)

Similarly we can derive G s,z from equations 2.50 and 2.52, leading to:

G s,z = 2 × 2 1 -ω ε e 1 + e 2 2 + 1 2 2 + 2α 2 e 2 + e 1 (α 1 -1) 4 × 2 1 -ω z e 1 +e 2 2 4 + 2 + α 1 e 2 + (1 + 2α 2 )e 1 α 1 (1 + 2α 2 ) -1
(2.53)

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

Like FCC and BCC z cells, FCC z cell is anisotropic so the calculation of an additional G s,x y is necessary to account for it. Similarly to BCC z , due to the additional strut on the side of the cell, we slice it into seven different parts as shown in figure 2.5b. To be consistent with both G s,x y calculation for BCC z cell and Yang et al. [22] method, we will only account for the influence of the additional strut added in the out-of-plane direction and not for the other struts located in slices S 1 and S 7 (see figure 2.5b). Actually, Yang et al. did not account for the influence of the struts that are perpendicular to the heat flow as they are not supposed to be relevant in heat conduction [22]. However, we want to differentiate BCC from BCC z and FCC from FCC z so we will account for their difference: the additional strut in the out-of-plane direction. Even if this seems paradoxical, it has proven to be correct (see part 2.3.3). From it, we derive the expression of P s 0 1 A s (s) d s using:

r 0 1 A s (s) d s = 2L/β A t ,1 + 2L -t 1 2r π , ( 2.54
) Using the definition of both G s,z and G s,x y for all cells and substituting them into equation 2.7, we can calculate for all of them k z and k x y and defined the so-called thermal conductivity tensor K e f f . The latter can be used in simulations to model the lattice structure infused with PCM and account for their anisotropy,

t 1 /2 r 1 A s (s) d s = t 1 -2 2 β L 2A t ,1 , (2.55)
(P s -t )/2 t /2 1 A s (s) d s = L -ω z t 2A l , (2.56)
K eff =     k x y 0 0 0 k x y 0 0 0 k z     .
(2.59)

Dimensionless parameters and constants determination

In order to calculate and compare the value of both G s,z and G s,x y we need to determine the value of all dimensionless parameters (i.e α i and e i ) and calibration constants (i.e As mentioned in chapter 1, determining the node size is of crucial importance to precisely estimate both porosity and effective thermal conductivity. The issue has been largely discussed, especially by Dai et al. [25] and Yang et al. [23]. The latter insisted on the geometry impossibility linked to a biased estimation of term e in both Boomsma and Poulikakos [24] and Dai [25] models. Actually, Dai corrected Boomsma and Poulikakos model to better fit experimental data in part by changing e value but pointed out that both models could not avoid geometrical impossible results. Yang indicated that using Boomsma and Poulikakos or Dai expression of porosity for foam would lead to the following statement : "as the porosity decreases, the ligament radius would increase while the length of the node remains constant to a point where the former exceeds the latter leading to impossible geometrical results". To improve the model Yang used a third degree polynomials [23] which parameters were determined to better fit experimental results from Calmidi and Mahajan [26]. The expression of that last parameter is defined in equation 1.95.

In the case of lattice structures, due to their regular and periodic geometry a more systematic method is used. Vaissier et al. [10] studied several lattice structures to be used as heat exchanger. To efficiently calculate both solid volume fraction and specific area (area per volume unit) Vaissier calculated the node volume and area of a variety of lattices and expressed it as a function of the strut diameter, which is known. We used Vaissier's work to determine the nodes volume of the structures depicted in figure 2.1 and derive the value of t and A t and infer e i and α i . Vaissier determined that the node volume is proportional to the cubic strut diameter and defined a proportionality ratio (γ). For a given node we can write: • BCC node (see figure 2.6a): This type of nodes is found on the corner and the center of a BCC cell, on the corner of a FCC cell (index 1) and at the center of a BCC z cell (index 2).

A t t = γ(2r ) 3 . ( 2
• BCC z node (see figure 2.6b): It is found on the corner of both BCC z (index 1) and FCC z (index 1) cells.

• FCC node (see figure 2.6c): This last node type is found at the center of both FCC (index 2) and FCC z (index 2) cell sides.

Due to their respective shape, we can divide the three node types into two categories: cube and cuboids. Actually, while both BCC and BCC z nodes dimensions are relatively similar along all three directions of space, FCC node dimensions along x and z (≈ 2(2r ))

axis are superior to the one along the y axis (≈ 2r ). Following this approximation, we represent both BCC and BCC z nodes as cubes of side t with A t = t 2 . FCC nodes are modeled as cuboids of height t and having a base of area A t = 2r t . We remind that A l = πr 2 and rearrange equation 2.60 to draw the expression of α:

α = 4 π γ 2/3 , (2.61)
for cubic nodes and

α = 4 π γ 1/2 , ( 2.62)
for the unique cuboid one.

Dealing with the expression of e we need to differentiate two cases: starting from equation 2.60 we see that the product A t t is expressed as a function of r . To link e = t /L to r we need to find a relation between r and L. Yet, we already defined β = P s /r and L can be related to P s through basic geometry relations that depends on the cell we consider: for a cuboid node in a FCC or FCC z cell. We notice that no matter the expression of e it always depends on β which is directly linked to the porosity. This agrees with Yang et al.

L = 3/
[23] conclusions stating that a constant value of e would lead to geometrical impossible results. We report all dimensionless parameters expressions and values in table 2.1 for each cell for more clarity.

Before G s,z and G s,x y can be calculated, the different calibration constants (i.e ω ε , ω z and ω x y ) need to be evaluated. For that, we first run a series of simulations on the four different cells presented here. Those simulations aim at determining the ETC for each cell at different porosity levels. The value of the calibration constants will be set for the model to fit the simulation data.

The above mentioned simulation can be described as follows: five cells are stacked up.

On top and at the bottom of those ones two plates, called heat spreaders, are placed for heat flow uniformity reasons. An example of such structure is given in figure 2.7 for BCC cells. Before actually running the simulation using COMSOL MULTIPHYSICS™, we first measure the actual porosity for each cell type to calibrate ω ε . Then, we run the simulation and apply a given heat flux (q) on top of that structure and measure the mean temperature on the top and on the bottom surface of the central cell. The choice of this particular set of temperature measurements will be detailed in part 2. temperature (∆T cc ) and express the ETC of the considered cell using Fourier's law:

k cc = q P s ∆T cc . ( 2.66)
Using this method, we are able to theoretically measure the ETC on a given porosity range for each cell and use those results to determine the calibration constants. Results

of the simulation and model fitting are displayed on figure 2.8. We report the values of the different calibration constants in table 2.1. On figure 2.8 we represent the normalized ETC with respect to the solid thermal conductivity (k s ). We notice an almost perfect fit between the analytical model and simulation data for porosity ranging from 0.55 to 1 which theoretically validates the model. This will be further discussed in the next part but we already notice that for a given porosity FCC z has the highest ETC in the out-of-plane direction while in the in-plane direction BCC behaves better than the other cells. This is mainly due to strut orientation and strut radius for a given porosity. It is pretty obvious that orienting struts in a given direction will help conduct heat in that particular direction explaining the great behavior of the FCC z cell in the out-of-plane direction. Actually, the struts on the cell sides are oriented in the out-of-plane direction and the other form a 45

•
angle with it which is lower than the 55 • angle from most struts in the BCC or BCC z cells. On the other hand depending on strut arrangement, the strut radius for a given porosity will change. This fact is easily illustrated when comparing, for example, a FCC or a FCC z cell:

adding a strut to the FCC cell to obtain a FCC z one at constant strut radius will obviously decrease the porosity as more conductive material is found in the CV. As high cross section is key to a good heat conduction adding more strut might not be profitable in terms of ETC. To illustrate that fact, we depict in figure 2.9 the corresponding porosity for a given given pore size. We observe that FCC is the one that provides the highest porosity level for a given strut radius which can be interesting both in term of ETC improvement but also to overcome additive manufacturing limitations.

Discussion

To asses ETC improvements from foam to architected structures, we need to find a way to properly compare them. Given the number of experimental data that can be found in the literature [27; 28] and the variety of different materials (copper, aluminum, plastics, nickel...) and manufacturing methods used to produce the foam it might not be relevant to compare their ETC or even their normalized ETC as the fluid used to fill the porous material might differ. Instead of directly comparing the ETC of the different structures for which we adapted Yang's model, to the foam, it would be more relevant to compare both G s,z and G s,x y respectively. Actually most of the ETC depends on the solid structure topology and material. We compare G s,i parameters of the lattices to the one of the foam and we add data extracted from experimental measurement from Ranut [27], who performed an extensive review on foam ETC modeling and measurement as well as Yao [28] who performed himself measurements on copper foams, to illustrate Yang's model preci- Aluminum foams [27] Copper foams [28] (a) Aluminum foams [27] Copper foams [28] (b) in figure 2.10a we observe an up-to 75% increase in G s,z from foam to FCC z , up-to 53% increase from foam to FCC, up-to 30% from foam to BCC z and up-to 7% increase from foam to BCC. In the mean time we observe, in figure 2.10b, an up to 7% increase in G s,x y from foam to BCC (BCC cell is isotropic) while almost no improvement nor decline in G s,x y from foam to BCC z . In opposition to the BCC based cells, FCC shows a decrease of about 13% and FCC z a decline of 25% in G s,x y . Both increases and decreases are linked to strut orientation and position both acting on the reduction or increase of the thermal path that the heat flow has to follow to go from one side of the cell to the other. It should be noted that while Yang's model validity requires the porosity to be over 0.9, model presented in this paper has been proven to be valid for porosity as low as 0.55 through simulations (see figure 2.8). This allows more freedom in the design of the thermal management device to better adapt to tougher working conditions, especially high power density components requiring higher ETC to cope with higher heat flux. Finally, we can also add that when considering structure manufacturing we need to take into account the limitations of additive manufacturing. Most metal additive printer use a powder bed locally heated by a laser spot that will melt the powder grain at the desired locations. Both powder grain diameter, laser spot power and diameter will affect the precision but at some point may also affect mechanical integrity. For these reason and using aluminum (AS7 alloy) strut diameter limit is set at 500 µm by all manufacturers that have been contacted.

Considering both G s,z and G s,x y value but also manufacturing limitation FCC cell seems to be the optimum one. It actually combines the second highest G s,z and the third G s,x y both by small margins. But the main advantage of the FCC cell is high strut diameter for a given porosity. This will allow more freedom in terms of porosity range allowing a better adaptability to a variety of problems. This is particularly important as in some cases a relatively tight space is available limiting the number of cells that can be stacked up making homogenization difficult as side effects that are usually neglected in the case of foam cannot be anymore.

We may now have to consider side effects that are not taken into account by this model in which we consider an isolated cell. Seeing stacked up cells and their environment as a whole as it is done on figure 2.7 we observe a brutal change of cross section between the so-heat spreader and the structure struts leading to thermal constrictions [29] (heat flow tightening) which limits heat slow may reduce effective thermal conductivity. This phenomenon is described in the next part.

Constriction Phenomena

Constriction General Definition and Equations

For mechanical purposes and to avoid any PCM leaking, the composite material formed by the PCM and an array of the selected cell need to be packaged. This is achieved by simply printing walls around the composite. Those walls have two purposes: avoid PCM from leaking and serve as thermal interface or heat spreader between the electronic com- that should be similar to the one defined on equation 2.66 as the exact same structure is modeled. The ETC (called k hs ) is calculated using equation 2.67 using data from a series of simulations in which pore size ranges from 1 to 10 mm, both heat spreader thickness is fixed at 1 mm while porosity is fixed at 0.8:

k hs = q 5P s ∆T ht . (2.67)
Unlike what is predicted by equation 2.7 we notice, in figure 2.12, that the ETC varies with pore size yet the porosity and topology is fixed, when considering the set of temperatures close to the heat spreaders while it is not when considering the set of temperatures close to the central cell. This is the sign that side effects are changing heat flow behavior at the interface and needs to be modeled for a better precision in the homogenization process.

This phenomenon is used to explain a variety of interfacial issues at the macroscopic scale (non-perfect contact between flat or rough surfaces (see figure 2.13), perfect contact between solids of different cross sectional area) or at the microscopic scale. It has

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

been described in a variety of fields including electrical engineering [30; 31] where the real area of contact was described. Actually, the interface between two solids in contacts is usually far from being flat and resemble more what is described in figure 2.13 : a series of micro-contacts created by the two solids roughness. Such interface is described by Greenwood [31] as being composed of number of clusters of micro-contacts which position is determined by the large-scale waviness of the surface and the micro-contacts by the small-scale surface roughness. Around those micro-contacts we notice a wide reduction of the cross sectional area leading to an additional electrical (and we will see later thermal)

resistance. Due to the similarity of both electrical and thermal phenomena, constriction is also observable for the heat flux at the interface between two solids. It is well depicted in figure 2.13 as when the heat flux get to the interface, it tends to avoid the low conductive area (vacuum on this figure) to go through the micro-contacts. The heat flux line bending creates the additional thermal resistance. This phenomenon occurs at many scale and while most models were developed to described what is happening at the micro-scale, some also deal with the issue of clusters [31; 32] or the perfect contact between solids of different cross sectional area [33][34][35] at the macro-scale or even at the nano-scale to describe the influence of phonon scattering between flakes or grains [36][37][38]. In a similar way Degiovanni [39; 40] studied the combination of several scales of roughness and their influence on constriction. The influence of the shape (triangular, circular, square...) [41][START_REF] Yovanovich | Spreading resistance of isoflux rectangles and strips on compound flux channels[END_REF][43] of both micro contact and heat spreader is also dealt with as this obviously influence the way heat will be conducted and so influence the value of the thermal constriction resistance.

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

Event though, as we have seen, the shape on both micro-contact and heat spreader have an influence on the value of the constriction resistance, a pretty convenient representation of the thermal contact could be represented as a cylinder of radius b and height h on which a heat flux q is applied on a micro-contact represented by a circle of radius a (see figure 2.14). Depending on the boundary conditions applied to each of the surfaces of the above mentioned cylinder, heat distribution and therefore thermal constriction resistance could be highly different. Gladwell [29] proposed an extensive review on that subject. To derive the value of the thermal constriction resistance, Gladwell considered the above mentioned cylinder to be surrounded by a fluid at temperature T ref on which is applied a total heat flux:

Q = 2π a 0 r q(r, 0)d r. (2.68)
To determine the expression of the constriction thermal resistance, Gladwell compared the difference between the mean temperature over the contact spot on top of the cylinder denoted S 1 (see equation 2.73) if a = b or if a < b. We note that temperature difference T c . Using it we can define the general expression of the constriction resistance

R c : R c = T c Q . ( 2.69)
For practical reasons Gladwell also defined a dimensionless constriction resistance Ψ, normalized with respect to the value 1/(4ak s ) related to the isothermal contact spot on the otherwise insulated surface of a half-space (see case a.2. below):

Ψ = 4k s aR c . (2.70)
Let Ω denotes the interior volume of the cylinder defined in figure 2.14. To determine the expression of T c we solve the following equation on Ω:

∇ 2 T(r, z) = 0, on Ω, (2.71)
Figure 2.13: Non-perfect micro-contact between two solids due to their respective roughness [44]. with the following general boundary conditions:

k s ∂T ∂z -H 1 T = f 1 (r ), on S 1 , k s ∂T ∂z -H 2 T = f 2 (r ), on S 2 , k s ∂T ∂r + H 3 T = f 3 (z), on S 3 , k s ∂T ∂z + H 4 T = f 4 (r ), on S 4 , (2.72)
where the (H i ) i ∈ [1,4] denotes convective exchange coefficient while the ( f i ) i ∈ [1,4] are real functions such that T(r, z) is continuous, and T(r, z) → 0 as r 2 + z 2 → ∞. They are defined on their respective surface of the cylinder (S i ) i ∈ [1,4] expressed as:

S 1 = {(r, z) : 0 ≤ r ≤ a, z = 0}, S 2 = {(r, z) : a < r ≤ b, z = 0}, S 3 = {(r, z) : r = b, 0 ≤ z ≤ h}, S 4 = {(r, z) : 0 ≤ r ≤ b, z = h}, S 0 = S 1 ∪ S 2 .
(2.73)

As pointed out by Gladwell, equation 2.71 could be solved using the general boundary conditions defined in equation 2.72 by superimposing basic solutions. Instead of that, Gladwell describes a series of particular cases that may represent specific situations encountered when dealing with micro contacts and will allow simplifying calculations. Six situations are categorized according to increasing complexity. We quickly described the cases defined by Gladwell, which mathematical definition is available in the literature, only detailing cases a and d as the former serves as the reference to normalize the thermal constriction resistance while the latter will be implemented in the ETC model.

a) The radius b and the height h of the cylinder are large compare to the contact spot radius a. The system is represented as a half space in which S 2 is insulated. and Yovanovich [44]. In this case it yields:

q = Q πa 2 , r < a; T(r, 0) = 2Q π 2 ak s E r a , r < a; Ψ = 32/3π 2 = 1.0808, (2.74)
with E(r /a) an elliptic integral of the second kind.

2) Instead of specifying the heat flux on S 1 the temperature is specified and set at T 1 . This time, it is considered a MBVP increasing the complexity of the solution.

The solution itself is due to Titchmarsh [45] and Sneddon [46]. In this case:

q = (2k s T 1 /π)(a 2 -r 2 ) -1 2 , r < a; T = T 1 , r < a; Ψ = 1. (2.75)
Those cases will serve as reference to normalize the dimensionless thermal constriction resistance due to the relative simple geometrical definition.

b) In this second case the cylinder is still a half-space but boundary conditions on S 1 and S 2 are two members of the set:

T = T 0 , ∂T/∂z = 0, k s ∂T.∂z = H(T -T 1 ). (2.76)
Situations that were analysed by Gladwell et al. [47] and for which Lemczyk and

Yovanovich [48] gave numerical results. These problems can be described as specific cases of the next ones described.

c) In the next case the cylinder radius b is finite while its height h is infinite. We consider S 2 to be insulated while S 3 could be either either insulated or isothermal. Once again two cases are considered:

1) We consider a first case where the flux is fixed on S 1 . In this situation the heat flux is fixed on the entire surface S 0 leading to a relatively simple problem as the boundary conditions are unmixed on S 0 . This problem was solved by Gladwell [29] and involved the computation of an infinite integral.

2) A MBVP arises when instead of fixing the flux on S 1 , the temperature is fixed.

This problem was solved by Cooke and Tranter [49] who considered isothermal conditions on S 3 . They solved the problem using Dini series yielding:

T(r, z) = ∞ n=1 σ -1 n a n J 0 (σ n r ) exp(-σ n z), ( 2.77)
where the σ n are the positive zero of J 0 (σ n b) = 0, J 0 being the Bessel function of order zero and of the first kind.

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

d) Gladwell [29] then described another unmixed problem in which all dimensions of the cylinder (b, h) are finite. Regarding the boundary conditions: S 3 is insulated while the heat flux is fixed on S 1 . S 3 and S 4 could either be insulated, isothermal or subject to general convective conditions. This case is detailed below as it corresponds to the problem encountered at the interface between the lattice structures and the heat spreaders.

e) Similarly b and h are considered finite. The region S 1 is subject to a convective condition, S 2 is either insulated or isothermal while S 3 and S 4 have general convective boundary conditions. This case is dealt with by Gladwell [29] but will not be detailed here.

f) In the last case, the dimensions of the cylinder are once again finite while all surfaces are subject to general convective boundary conditions. This problem was also solved by Gladwell [29] and will not be detailed here.

We choose here to detail case d as this is the one that resembles the more the configuration that needs to be modeled in figure 2.7. Actually, the structure represented in figure 2.7 could be considered as a pattern of a network spread across the plan defined by the heat spreader. If the network spreads across a sufficiently large surface, the study of a single pattern can give a proper approximation of the behavior of the entire network.

Profiting from the symmetry of the structure, we consider that its sides are thermally insulated while the heat flux is fixed on the top surface and the bottom surface is subjected to general convective boundary conditions. Concentrating on one of the two heat spreaders, starting with the top one, we actually notice that a given heat flux enters the top of it and if the PCM thermal conductivity is low enough compared to the solid one, it can considered that the entire heat flux has to come out of it through a spot which surface corresponds to the one of a node of the cell. This situation is close to case d and it will serve as a reference.

For that reason we detail the way it is solved.

To do so we need to re-write the boundary conditions on the different surfaces S i of the cylinder originally defined in equation 2.72 as follows:

k s ∂T ∂z = Q πa 2 , on S 1 , ∂T ∂z = 0, on S 2 ,
∂T ∂r = 0, on S 3 ,

∂T ∂z = A 1 , on S 4 , (2.78)
in which A 1 is defined by the equation:

Q + πb 2 k s A 1 = 0, (2.79) CHAPTER 2.

PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

stating that all the heat entering Ω through the contact spot of radius a must leave it through the bottom surface.

A general solution of equation 2.71 was given by Sneddon [50] yielding: To determine the integration constant we re-write the partial derivative of T at the required locations (i.e. locations defined in equation 2.78):

T(r, z) = A 0 + A 1 z + ∞ 0 Γ(ξ) cosh(ξ(h -z))J 0 (ξr )d ξ + ∞ 0 Λ(ξ) sinh(ξ(h -z))J 0 (ξr )d ξ + ∞ n=1 λ -1 n C n cos(λ n z)I 0 (λ n r ) + ∞ n=1 λ -1 n D n sin(λ n z)I 0 (λ n r ), ( 2
∂T ∂z r,0 = A 1 - ∞ 0 ξΓ(ξ) sinh(ξh)J 0 (ξr )d ξ - ∞ 0 ξΛ(ξ) cosh(ξh)J 0 (ξr )d ξ + ∞ n=1 D n I 0 (λ n r ), (2.81) ∂T ∂z r,h = A 1 - ∞ 0 ξΛ(ξ)J 0 (ξr )d ξ - ∞ n=1 C n sin(λ n h)I 0 (λ n r ) + ∞ n=1 D n cos(λ n h)I 0 (λ n r ), (2.82) ∂T ∂r b,z = - ∞ 0 ξΓ(ξ) cosh(ξ(h -z))J 1 (ξb)d ξ - ∞ 0 ξΛ(ξ) sinh(ξ(h -z))J 1 (ξb)d ξ + ∞ n=1 C n cos(λ n z)I 1 (λ n b) + ∞ n=1 D n sin(λ n z)I 0 (λ n b), (2.83)
using the following identities :

∂ ∂r J 0 (ξr ) = -ξJ 1 (ξr ), ∂ ∂r I 0 (ξr ) = ξI 1 (ξr ) (2.84)
To satisfy the boundary conditions on S 4 , equation 2.82 requires Λ(ξ) = 0, D n = 0 and si n(λ n h) = 0. We derive:

λ n = nπ/h, ∀n ∈ N + . (2.85) CHAPTER 2.

PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

Now, using the boundary conditions on S 0 , we have:

∞ 0 ξΓ(ξ) sinh(ξh)J 0 (ξr )d ξ = A 1 + q(r )/k s , if r < a A 1 , if a < r < b. (2.86)
If we write Γ(ξ) = Γ 1 (ξ) + Γ 2 (ξ), equation 2.86 requires:

Γ 1 (ξ)si n(ξh) = A 1 bξ -1 J 1 (ξb), (2.87)
We actually note that, in this case:

∞ 0 ξΓ 1 (ξ) sinh(ξh)J 0 (ξr )d ξ = A 1 ∞ 0 bJ 1 (ξb)J 0 (ξr )d ξ = A 1 (2.88)
Using the following identity [51] with a = r , b = b, µ = 0, ν = 1 and λ = 0:

∞ 0 J µ (at )J ν (bt ) t λ d t = a µ Γ 1 2 (ν + µ -λ -1) 2 λ b µ-λ+1 Γ 1 2 (ν -µ + λ + 1) × F 1 2 (µ + ν -λ + 1); 1 2 (µ -ν -λ + 1); µ + 1; a 2 b 2 . (2.89)
Finally the expression of Γ 2 is given by:

∞ 0 ξΓ 2 (ξ) sinh(ξh)J 0 (ξr )d ξ = q(r )/k s , if r < a 0, if a < r < b. (2.90)
Using the same identity we can find:

Γ 2 (ξ)si nh(ξh) = Q πak s J 1 (ξa) ξ . ( 2

.91)

We now use the insulated boundary condition on S 3 which leads to:

∞ n=1 C n I 1 (λ n b)cos(λ n z) = ∞ 0 ξΛ(ξ)J 0 (ξr )d ξ. (2.92)
To derive the value of C n we first use the orthogonality of cosh(λ n z) over [0, h], leading to:

C n I 1 (λ n b) = 2 h ∞ 0 ξΓ(ξ)J 1 (ξb) h 0 cos(ξ(h -z))cos(λ n z)d zd ξ. ( 2

.93)

with:

h 0 cos(ξ(h -z))cos(λ n z)d z = ξsi nh(ξh) λ 2 n + ξ 2 .
(2.94)

Then we have:

C n I 1 (λ n b) = 2 h ∞ 0 ξ 2 si nh(ξh) λ 2 n + ξ 2 Γ(ξ)J 1 (ξb)d ξ.
(2.95)

To calculate the integral on the right side of the equation, we substitute both equations 2.87 and 2.91 in equation 2.95 in combination with the following identity [51] in which K 1 is the modified Bessel function of the second kind:

∞ 0 ξJ 1 (ξa)J 1 (ξb) λ 2 n + ξ 2 d ξ = I 1 (λ n a)K 1 (λ n b).
(2.96)

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

Finally we derive:

C n I 1 (λ n b) = 2 h A 1 bI 1 (λ n b) + Q πak s I 1 (λ n a) K 1 (λ n b).
(2.97)

We are now able to fully express the temperature over the contact spot. We note that when a = b the temperature over the entire S 0 surface can be expressed as T(r, h) = A 0 . We can now express T c as:

T c = 2 a ∞ 0 ξ -1 Γ(ξ) cosh(ξh)J 0 (ξa)d ξ + ∞ n=1 λ -2 n C n I(λ n a) . (2.98)
We note:

cosh(ξh) = 1 + ν(ξh), (2.99)
with:

ν(s) = 2e -2s 1 -e -2s .
(2.100)

Using the expressions of Γ, and C n we are now able to write the expression of the dimensionless constriction resistance in equation 2.70 as a sum of three different terms:

Ψ = θ + φ + χ, (2.101)
with:

θ = 8 πa ∞ 0 ξ -2 J 1 (ξa) J 1 (ξa) - a b J 1 (ξb) d ξ, (2.102) φ = 8 πa ∞ 0 ξ -2 ν(ξh)J 1 (ξa) J 1 (ξa) - a b J 1 (ξb) d ξ, (2.103)
χ = 16 πha ∞ n=1 λ -2 I 1 (λ n a) K 1 (λ n b) I 1 (λ n b) I a (λ n a) - a b I 1 (λ n b) . (2.104)
In equation 2.101, θ only depends on β = b/a and could be seen as an image of how much the heat flux lines have to bend to verify the boundary conditions. We actually see in figure 2.15 that as a → b this term tends to vanish. The second term φ depends on both β and η = h/a: it is an image of how much space the heat flux lines are given to bend.

It seems that increasing this space actually reduces the constriction resistance as φ → 0 when h → ∞. The last term χ depends on the boundary conditions on S 3 a S 4 . It is an image of the heat loss that could occur on S 3 depending on the boundary conditions on this surface. In the case where S 3 is insulated χ is negligible compare to φ and θ and will be neglected in the following parts.

Implementation in the Effective Thermal Conductivity Model

In the previous section we have been able to determine the expression of the thermal constriction resistance created when heat flow through a relatively small contact spot into a

θ = 8 π α 1 r ∞ 0 ξ -2 J 1 (ξ α 1 r ) J 1 (ξ α 1 r ) - α 1 A s /π P s / π J 1 (ξP s / π) d ξ, (2.105)
φ = 8 π α 1 r ∞ 0 ξ -2 ν(ξh)J 1 (ξ α 1 r ) J 1 (ξ α 1 r ) - α 1 r P s / π J 1 (ξP s / π) d ξ. (2.106)
To explain the phenomena observed in section 2.3.1, we need to implement thermal constriction resistance in the ETC model. To do so, a structure similar to the one displayed in figure 2.7 but involving N cells stacked on top of each other instead of five is represented. Using the assumption made in section 2.2.1 of parallel heat flow between the PCM and the porous material, we represent the above mentioned structure as a network of thermal resistance depicted in figure 2.17 with:

R s = P s A 0 k s G s,z (1 -ε) , (2.107)
and

R f = P s A 0 k f ε . (2.108)
Using this representation we are able to derive the ETC of such structure including side effects due to thermal constriction:

k z = k s (1 -ε)G s,z 1 + 2R c A 0 G s,z (1 -ε)k s P s N -1 + k f ε. (2.109)
This new expression of the ETC helps explaining the apparent dependency of the ETC on the pore size which did not appear in the previous models especially those considering foams. It appears that this dependency is only caused by side effects that are seen when only few cells are stacked on top of each other. This could also explain why most studies on foams that searched the influence of the pore size showed a small influence of this parameter. The influence of the side effect declines quite fast (∝ 1/N) when the number of

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

cell increase, yet most models calculating foam ETC are made to account for the statistical behavior of such structure which implies considering a high number of cells making side effects negligible.

Numerical validation and Discussion

To justify the use of ETC model defined in equation 2. [START_REF] Mahdi | Solidification enhancement with multiple PCMs , cascaded metal foam and nanoparticles in the shell-and-tube energy storage system[END_REF] for constriction seem to improve the prediction. We observe that as the pore size decreases and draws near the heat spreader value, ETC increases. This fact is explained by the evolution of β and η in this situation: if we remind that P s = βr we notice that when P s draws near h the same is true for the strut radius which tends to decrease the value of η hence the one of the constriction resistance. Nonetheless, it is important to note that ETC dependency of pore size is uniquely due to the presence of constriction which means that if the structure considered in the simulation was made of a sufficiently high number of stacked-up cells or if we considered an isolated cell none of this effects would be observable. This model aims at making easier design and manufacturing of PCM architected thermal conductivity enhancers. It helps determining optimum cell topology as well as best surroundings to avoid a too strong influence of the side effects. Model depicted on equation 2.109 shows dependency of the ETC on five topological parameters that can be tuned using additive manufacturing:

R c R c

NR

• Cell general topology (BCC, FCC...),

• Porosity (ε),

• Pore size (P s ),

• Heat spreader thickness (h),

• Number of stacked cells (N). First, as FCC and FCC z provide higher ETC it seems that orientation of the the struts closer to the heat flow direction is beneficial to the ETC. Actually, for both BCC and BCC z structures, most struts have a higher angle from the heat flow direction, increasing the heat path from one side to the other side of the cell. It is interesting to observe that straight fins oriented in the heat flow direction would provide G s,z = 1: in that case the law of mixture would be valid in that particular direction leading to a maximum ETC but the inplane ETC would be drastically reduced. On the other hand, one can observe that for all cell topologies ETC increased with a decreasing porosity in a non-linear way unlike it is planned by law of mixture usually used to determined composite properties [27]. It was also proven that for a given porosity and heat spreader thickness, increasing pore size would be detrimental to the ETC (see figure 2.18). The same outcome is observe when decreasing the heat spreader thickness for a given porosity and pore size (see figure 2. 19). Both statements are actually linked to the same phenomenon: thermal constriction is more dependent on those dimensions ratio than on their actual value which is clearly seen in equations 2.102 and 2.103. Actually, even if it is less practical in terms of manufacturing, describing structure by the value of β = P s /r and η = h/r would give more information on their ability to conduct heat. This is visible in figure 2. 19 where the dimensionless ETC is expressed as a function of the heat spreader thickness for 5 values of strut radius with fixed value of β. Here we observe that as long as the heat spreader thickness is superior to r α 1 (which corresponds to the node radius that is in contact with the heat spreader), the heat spreader thickness has almost no influence on the ETC.

Finally, when the number of cells in the heat flow direction increases, the value of both constriction resistances on top and at the bottom of the cell 2 × R c tends to become negligible compared to the total thermal resistance of the porous structure N × R s as seen in figure 2.17. In that particular case, both models in equations 2.7 and 2.109 are equivalent.

Experimental Validation

Sample description and manufacturing

To further deepen the study of the ETC of architected structures we want to experimentally validate the models in equations 2.7 and 2.109 by measuring the effective thermal ameter which allowed porosity to range from 0.73 to 0.9 depending on the cell type and pore size using Selective Laser Melting (SLM). The SLM process consists in building a given structure from a Computer Aided Design (CAD) file by locally melting parts of powder layers thanks to a laser beam [52]. The building process starts with a powder layer on a substrate that is melted by a high energy density laser beam at selected location defined in the CAD file. Once the first layer has been melted at the right locations another layer is added on top of it and melted at selected locations again. The process is repeated until the entire structure is printed. After this steps, the un-melted powder is removed to obtain the desired design. Many parameters can be tuned to adapt to the searched resolution and powder material such as:

• Powder particle size : larger particles lead to a loss of resolution while smaller ones can agglomerate together through van der Walls forces resulting in poor powder flowability and poor powder deposition. In that sens powder particle size usually ranges from 20 to 100 µm,

• Laser wavelength which has to be wisely chosen to maximize photon absorption by the powder. Most SLM machines today use Yb:YAG fibre laser with a wavelength λ =1.06 µm [52],

• Volumetric energy: this energy is the one used to heat up and melt powder particles and is characterized by three parameters that can be tuned : laser power, scanning speed and layer thickness (see figure 2.21). The combination of these three parameters will affect molten pool wetting with the previous layer, provoke material evaporation, cause keyhole effect (which is due to vapor forming in the molten powder causing porosity [53]), induced porosity in the printed material.

For this series of sample, LISI-AEROSPACE™ used 50 µm layers with a powder of AS7 aluminum (more precisely AlSi7Mg0.6 alloy). It gave decent results (see photo of all samples in table A.1 in appendix A.1.2), however some defects are observable. First, the surface is not perfect, as expected, and roughness is observable on some samples (see fig-

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

• width: 5.47 mm,

• height: 5.61 mm.

Then, its mass is measured to be 5.61 g. From those measurements we derive the material density at room temperature (i.e. 20 • C) ρ = 2542 kg•m -3 given with a 10% uncertainty.

We, then, go on with the specific heat capacity measurement which is performed using a calorimeter (Calorimètre Calvet C80II from Setaram). Compared to other method like Differential Scanning Calorimetry (DSC), using a calorimeter allows the use of a heavier sample which ensures material uniformity. In addition, this method involves isothermal steps rather than dynamic ones which prevent any influence of the heating rate on the measurement [54]. This machine is composed of two identical crucibles : one is called the reference crucible and stay empty for all measurement while the sample to be measured is put in the second crucible. A series of thermocouples are used to determine the heat flux from the machine to both crucibles. Two sets of measurements are then performed: one with both crucibles empty (the signal obtained here would be withdrawn from the signal obtained in the second set of measurements) and the other one with a sample of mass m s .

For each set of measurements and for the different temperatures at which the specific heat capacity needs to be determined a temperature increase ∆T with a heating speed β is applied to both crucibles. We denote Q b and Q s the respective heat flux measured respectively in the first and in the second set of measurements from which we derive the expression of the specific heat capacity: The measurement itself is preformed as follows: both crucibles temperature are brought to the measurement temperature minus one degree Celsius. Following that, the crucible undergo a relatively long (here around 20 000 s) isothermal step which allows the temperature to be uniform in the entire sample. Then a temperature ramp of 2 • C with a speed of 0.05 • C•min -1 is performed. This step allows a relatively slow temperature rise for the temperature to be as uniform as possible in the sample to avoid any temperature gradient that could distort the measurement. Once done, another isothermal step is imposed to both crucible. During those different steps, both temperature and heat flux are constantly monitored to obtain a so-called enthalpic graph as shown in figure 2.23. From it, we derive the specific heat capacity as it corresponds to the air under the blue curve, ignoring a multiplication factor.

C p = Q s -Q b m s ∆T . ( 2
Before the measurement can be done on the actual sample we want to characterize, the calorimeter is calibrated using a POCO sample. POCO is a type a high quality graphite that is used for measurement machine calibration and electrode manufacturing due to its low particle size. This material has been characterized by a number of industrial and laboratory hence it can be used to calibrate the calorimeter. The results can be found in figure A.2 in Appendix A.1.3. On this graph the reference specific heat capacity values from the literature are within the 5% uncertainty margin compare to the one determined experimentally hence validating the calorimeter calibration.

Once the machine is calibrated, the measurement can be performed on the actual sample for five different temperatures : 25 Another parameter that is not directly related to the thermal diffusivity is measured: the coefficient of thermal expansion α. Even though, it is not related to the thermal diffusivity, the length of the sample plays a role in the derivation of the thermal diffusivity from the results of the flash laser method. To do so, Influtherm used a so-called dilatome- ter, which consists in heating a slab of the material to be characterized while precisely measuring the evolution of its overall temperature as well as the evolution of its length, the coefficient α is then obtained using the following expression:

α = d L L 0 1 ∆T , (2.111)
with:

• α: Coefficient of thermal expansion (K -1 ),

• L 0 : Slab length at room temperature (m),

• d L: Length increase at the considered temperature (m),

• ∆T: Temperature increase (K).

The elongation of the slab is measured using a TMA 402-F3 Hyperion that is able to work from -170 • C to 1000 • C. In this machine, the sample is held by a mobile part that apply a force on it, allowing it to stay in the right position for the measurement but that does not compress it so that the measurement is not skewed (see figure 2.24). When the slab is heated and starts to expand the movement of the mobile part is measured and the expansion of the slab can be derived. Several parameters can interfere in the measurement:

• Temperature gradient in the slab: as seen in figure 2.24, the thermocouple is roughly placed in the middle of the slab so this is the only spot where the temperature is measured. The slab dimensions are set to limit the gradient (as they are "small enough") but any temperature difference may induce mistakes in the measurement.

• Gravity: as we measure the expansion of the slab in the out-of-plane direction and because the expansion is small compared to the slab dimensions (α ≈ 10 -5 K -1 ) gravity may reduce the expansion of the slab as it might "fall under its own weight".

An alternative could be to measure the expansion in-plane but this would cause friction with the medium that may also reduce the expansion. A similar effect is induced by the mass of the mobile part.

• Heating rate: one must ensure that this parameter remains constant over the entire period of measurement as the coefficient of thermal expansion is derived from the evolution of the slab length and from the heating rate that is considered constant. Any change in the heating rate might induce error in the estimation of the coefficient of thermal expansion. This can be seen in figure 2.25 as at the end of the measurement the sample is subjected to an isotherm creating a virtual drop of the coefficient of thermal expansion for a temperature over 200 • C.

The coefficient of thermal expansion of the base material is represented in figure 2. 25 for temperature ranging from 20

Base Material Thermal Diffusivity measurement

The last parameter to be determined is the thermal conductivity. This parameter is determined indirectly using the flash laser method. Actually, flash laser method is a dynamic method that only gives an image of the thermal diffusivity. From its value and the ones of both density and specific heat capacity, the thermal conductivity can be derived. To do so, Influtherm used a machine of their own conception (see figure 2.26) that allows measurements under vacuum condition to avoid heat loss at temperatures ranging from -30

to 200 • C.
The thermal diffusivity of the base material is determined using a cylindrical sample with a height of 11.67 mm and a diameter of 24 mm. Its diffusivity is measured for five different temperatures (20, 50, 100, 150 and 190 • C) to determine if any sensitivity to this parameter exists.

Flash laser method consists in applying a Dirac of energy on the front face of the sample while measuring the temperature of the rear side. Using one of the many existing models in the literature [55][56][57][58], one is able to derive the thermal diffusivity of the studied sample for different experimental conditions that are accounted for through different boundary and initial conditions in the different models. Although any experimental conditions can be accounted for in the models, most of them include convective of radiative heat losses involving the estimation of convective or radiative exchange coefficient that are unknown for most materials. This motivates an experimental setup that avoid most heat losses. To do so, the sample is placed in a vacuum environment which reduces convective heat loss to the minimum, while radiative heat loss are reduced by the metallic enclosure the sample is placed in (see figure 2.26). This still leaves some ways for heat loss as the sample has to be placed on the setup, which creates heat loss by contact. Hence, the surface of contact between the sample and the setup is minimized as much as possible. Another challenge linked to the flash laser method lies in the temperature measurement. Actually, the different samples undergo small temperature variations (< 5 K) over a relatively short period of time. Any error on the temperature measurement would drastically reduce the precision of the method. Using an usual thermocouple to measure the lower face temperature would create a thermal contact resistance between the thermocouple and the sample itself causing innacuracy in the measurement. To avoid it, the solution used by Influterm involves a hot-junction-less thermocouple. As shown in figure 2.27, the hot-junction-less thermocouple virtually delete the thermal contact resistance as the sample itself plays the role of the hot junction. To do so, an electrical contact has to be maintained between the sample and the two metals forming the thermocouple which is done by pressing the sample on the two branches displayed in figure 2.26b.

The obtained thermogram (see figure 2.28) is processed using five different models (i.e. Parker [55], Clark and Taylor [56], Logarithm [59], Momentum [60], Partial Times [58]) in order to derive the thermal diffusivity of the sample. On the same thermogram, we compare the temperature profile obtained using the Parker method (Influtherm data) which displays an almost perfect fit and ensures that this method is suitable for this sample. From it, we derive the value of the thermal diffusivity for the five temperatures described earlier. We notice that the diffusivity remains constant with temperature. From it we derive the thermal diffusivity using the previously obtained specific heat capacity value while neglecting the variation of density with temperature. All value of thermal diffusivity and conductivity are reported in table 2.4. Those value are coherent with other values of thermal conductivity for similar aluminum alloy.

Those measurements (i.e. Density, Specific Heat Capacity, Coefficient of Thermal expansion) enabled to fully characterized the base material of the samples described in part 2.4.1 from a thermal point of view. These values can now be used to characterize the above mentioned samples which should validate the ETC models described in sections 2.2 and 2.3.

Adapted Flash Laser Method Results Processing

Now that the base materials of the samples described in table 2.2 is characterized, we need to characterize the ETC of those samples. To do so, flash laser method is used once again but before doing so, one should be reminded that flash laser method and the models aiming at deriving the thermal diffusivity are only valid for homogeneous, isotropic and solid samples.

Temperature Thermal Diffusivity Uncertainty Thermal Conductivity None of the above mentioned characteristic applied to the samples previously described. Actually, most samples are anisotropic and are, by nature heterogeneous. In addition, their radius to height ratio is highly detrimental to heat loss limitation which could cause heat conduction not to be unidirectional which is another of the key assumptions made in the different models of thermal diffusivity derivation. Although, it seems rather ill-advised to carry on flash laser on those samples, the advantages of the method (i.e. absence of thermal contact, dynamic method, heat loss limitation...) still makes it the most reliable method in this situation. Even though, the method can still be used, the assumptions made in the different models used to derive the thermal diffusivity are no longer verified which makes the above mentioned models unreliable in our case. does not falls on itself while being printed. This method is similar to the one used in the determination of contact resistance using flash laser method [61; 62]. We report the properties of the different layers as:

( • C) (m 2 •s -1 ) (m 2 •s -1 ) (W•K -1 •m -1 ) 25
• Heat spreader: k Al , ρ Al and C p ,

• Lattice structure:

k e f f , ρ e f f = (1 -ε)ρ Al and C p ,
• Interface:

k I = k e f f + k Al /3, ρ I = ρ Al (1 -ε + 1/3) and C p .
To determine the thermal conductivity of the interface, we consider a unit cell filled with two pyramids (at the top and at the bottom of it) representing a third of the cell volume. We then consider parallel heat flow between the struts and the pyramids. The mean cross section of the pyramid ( Ā) is defined, for a BCC structure, as: Using this representation of the samples, the following heat equation is solved by a finite difference scheme: -

Ā = 2 P s P s /2 0 (P s -2x) 2 d x = P 2 s 3 = A 0 3 . ( 2
1 α(x) ∂T ∂t x,t = ∂ 2 T ∂x 2 x,t - 2πh r πR 2 k(x) (T(x, t ) -T ext ) , ( 2
k(0) ∂T ∂x 0,t = Q(t ) -h 0 (T(0, t ) -T ext ) , (2.115)
with Q the heat flux due to the flash, it is time dependent as it is supposed to be short and h 0 is the convective exchange coefficient related to the top heat spreader. At the bottom of the sample, a simple convective exchange is considered leading to:

-k(L) ∂T ∂x L,t = h L (T(L, t ) -T ext ) . ( 2

.116)

Between each layer, thermal resistances are added to better account for the variations of average properties induced by the different topology from one layer to another. Two thermal conductances are defined: h b between the heat spreader and the interface and h t between the interface and the lattice structure. The value of h b is set relatively high compared to other conductance in the system. It actually represents the perfect contact between heat spreader and the base of the pyramids that are filling the pores close to it. Negligible thermal constriction is expected. On the other hand, h t is calculated using Cooper [63] model, who represented non-perfect contact of two solids by an array of cones in contact with each other. The expression was modified to account for the change of topology (from cone to pyramid). We define h t as:

h t = 4 1 k I + 1 k e f f -1 2.2 α 2 r P 2 s 1 - 2.2 α 2 r P 2 s -3/2 . ( 2

.117)

We now define the boundary conditions at the different interfaces separating each layer. We first define the B.C. between the top heat spreader and the first interface on the heat spreader side:

-k Al ∂T ∂x h -,t = h b T(h + , t ) -T(h -, t ) . ( 2

.118)

A similar B.C. is defined on the interface material side:

-k I ∂T ∂x h + ,t = h b T(h -, t ) -T(h + , t ) . (2.119) CHAPTER 2.

PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

The same set of equations is defined at the interface between the bottom heat spreader and the second interface zone. On the interface zone side the B.C. is defined as:

-k I ∂T ∂x (h+P s /2) -,t = h t T((h + P s /2) + , t ) -T((h + P s /2) -, t ) , (2.120)
while on the lattice structure side:

k e f f ∂T ∂x (h+P s /2) + ,t = h t T((h + P s /2) -, t ) -T((h + P s /2) + , t ) .

(2.121)

The model described here is combined with a gradient descent algorithm to determine the value of thermal diffusivity of the lattice and interface layers as well as the exchange coefficients that lead to the best fit of the experimental thermograms. Results for sample 3 and 5 are depicted in figures 2.29a and 2.29c and show great improvements compared to Parker model. For sample 3 residual has been greatly reduced while an almost perfect fit is obtained for sample 5. The difference of fitting quality can be explained by the structures of the two samples: samples 3 is composed of relatively large pores and have relatively high porosity compared to sample 5. This improvement in the experimental thermogram fitting provides a higher level of reliability in the thermal diffusivity and conductivity measurement, which was considered satisfying for all samples mentioned in table 2.2.

To precisely determine the ETC of the samples described in table 2.2, three of each sample type are manufactured to limit the influence of manufacturing on the measurements. In addition, the flash laser method is run for different temperatures to confirm the fact that thermal diffusivity is independent of temperature for porous samples. For each of these temperatures the flash laser method is run three to five times for repeatability purposes. We report the measurement conditions of each sample in table A.2 found in Appendix. In this table, we report the values of thermal diffusivity for each sample and each temperature with a 10% uncertainty obtained with the previously described 1D finite difference model. We notice that, no matter the temperature, the thermal diffusivity remains almost constant for a given sample type confirming the independence of this parameter regarding the temperature. Using the values of porosity (ε), thermal diffusivity (α), density (ρ) and specific heat capacity (C p ) measured in the previous paragraph, we derive the value of the ETC for each sample at room temperature using: We summarize the average values of ETC obtained for each sample type in table 2.5 for clarity. We remind here that, even for samples adapted to the flash laser method (i.e. isotropic, solid, homogeneous and with a high radius to height ratio) thermal diffusivity results are given with a 5% uncertainty according to Influtherm standard. For this modified version of the thermal diffusivity determination, we could not expect uncertainty below 10% and results have to be interpreted accordingly.

k e f f = αC p ρ(1 -ε). ( 2

Discussion and Model Validation

In table 2.5, experimental values from flash laser measurements are compared to the ones obtained using both models defined in equation 2.7 and 2.109. From it, we determine a maximum difference of 12.81% between models and experiment for sample 5 which tends to validate the use of both model for ETC prediction of PCM architected enhancers.

We now compare models on specific cases. Samples 1, 2 and 3 have been designed to see the influence of the pore size (4, 6 and 8 mm) on the ETC. Obviously, model defined in equation 2.7 predicts an identical ETC for the three samples as it does not account for constriction while the contrary is true for the constriction model. When looking at the experimental values, a change is observed and is similar to the one predicted by the model though it is smaller than expected from sample 2 to sample 3. This may be due to the presence of the pyramids that may reduce the effect of constriction increasing the space

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

given to the heat flux lines to be constricted. In this case the validity of the constriction model is verified.

Samples 4, 5 and 6 were designed to observe the influence of porosity. In that case both models behave the same way than experimental values though constriction-free model gives better results. This might be explained by the slight difference between the model and the actual topology: struts tends to be larger than expected, this propagates to nodes topology changing how heat is transfer inside the sample. This is especially true for larger struts (i.e. smaller porosity). This effect is also combined to the one of pyramids that may also reduce the effect of constriction. Both models seems to be valid in that case as the effect of constriction is not clear. We observe that the constriction-free model gives better results here even though both models are within the error margin (around 10%).

The third set of samples is composed of samples 7, 8 and 9. Those have been designed

to examine the consequence of a change of heat spreader thickness (0.5, 1 and 2 mm). As mentioned before, if the heat spreader thickness is superior or close to the strut radius its influence on the ETC is negligible. Due to both manufacturing limitations and measurements requirements on sample dimensions, heat spreader thickness is of the same order of magnitude as the strut radius limiting its influence. Difference predicted by the constriction model is small. As a consequence experimental and model values do not follow the same trend and it is hard to conclude on the model validity from an experimental point of view. However both constriction and constriction-free predict ETC values close to the ones obtained experimentally.

A final set of samples composed of sample 10, 11 and 12 has been designed to corroborate the influence of the number of cells in the heat flow direction on the ETC. For the same reasons than the previous set, the difference between each value is relatively small and may lie within the error margin. Even though experimental and model value fluctuate the same way it is hard to conclude on the validity on the constriction model from this experiment. However both models predict value close to the experimental ones.

Conclusion

In this chapter, the development and validation of a model aiming at predicting the ETC of architected lattice structures for the enhancement of PCM is detailed. Using this particular model, the efficiency of such structure relatively to the one of more usual enhancers (especially foam) is assessed and showed promising results. An up to 75% ETC gain was observed for some directions in space.

Based on the work of Yang et al. [22; 23], the mathematical formulation of the ETC model for architected structure is detailed on an isolated cell. In this part, a first set of

CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

simulations allowed a first theoretical validation of the model which showed good agreement. This first model was used to compare the different structures studied in this chapter to the most usual PCM enhancer: foam. Comparison showed that FCC z gave the highest ETC in the out-of-plane direction as a it is 75% higher than the one of the foam for a given porosity. On the other hand, BCC showed the highest in-plane ETC with a 7% increased compared to foam. This model and comparison show the potential of such structure in terms of PCM enhancement.

In a second part of this chapter, the impact of thermal constriction at the interface between the structure described in the previous chapter and a hypothetical PCM enclosure is assessed. Using a model developed by Gladwell [47], constriction was implemented into the above mentioned model and allowed a more refined ETC estimation. Using this improved model, additional parameters of the porous network are accounted for, like for example pore size, number of stacked cells or heat spreader dimensions.

The final part of this chapter aimed at experimentally validating the ETC model. To do so, a series of samples were designed to assess the actual influence of each of the previously cited parameters on the ETC of architected PCM enhancers. To measure the ETC of the different samples, the flash laser method was used. Although this method has shown great results in a number of situations, it is not adapted to porous, anisotropic and heterogeneous samples like the one we designed. In that sense, it was necessary to adapt the way flash laser thermograms were processed to derive reliable results. By developing a 1D finite difference model we were able to do so. Comparing the experimental results to the one predicted by both models showed good agreement as a maximum 12.81% difference was observed between models and experiments. In addition, while both models seemed to globally agree with the experimental results, the constriction model seemed to show more refined estimation in some cases.

In this chapter, we were able to develop and validate a model for the estimation of the ETC of architected structures that can be manufactured using 3D printing. Additive manufacturing (or 3D printing) offers great freedom in design which can be used to optimize topology as much as possible to enhance PCM based thermal management devices.

In that sens, the next chapter of this manuscript will concentrate on the use of concept called Cascaded PCM [64][65][66][67][68] which consists in adapting the filler ratio depending on the distance from the heat source to maximise heat transfer and limit temperature gradient to make the best out of the PCM heat absorption capacities.

Introduction

In part 1.4.1, through the study of Ahmed et al. [1] work, it was proved that the topology of the filler for a given porosity and material would influence the Effective Thermal Conductivity (ETC) of the composite. The question was further studied in chapter 2, where an ETC model for Lattice structure was developed and compared to the one existing for foams. It was shown that a change of topology could increase the ETC by 35% for a given porosity and material which was confirmed experimentally. To further confirm the interest of such filler, it is thought to be applied on a given use-case of industrial interest.

This use-case is constituted of four dies each generating 100 W for 30 s in an environment at 80 • C and which surface temperature must be maintained under 120 • C for them to remain functional. It should be added that due to the final end application, the Phase Change Material (PCM) discharge will not be studied here. To perform this, a cylindrical space with a radius of 60 mm and up to 11 mm is available. The current solution for such use-case is composed of a disk having the previous dimensions and made of aluminum alloy (Al6061). The current solution needs improvement in terms of thermal management as it cannot reach the 120 • C (see part 3.3.2) and in terms of mass for integration purposes.

Before the use-case can actually be studied, modeled and simulated, the optimum PCM to be used for such application needs to be determined and characterized. Actually, as seen in part 1.3.1, the enthalpy method can be used to simulate the fusion of the PCM but it requires a good representation of the evolution of the specific heat capacity as a function of the temperature. Usual method for the estimation of this parameter is the Differential Scanning Calorimetry (DSC), which was already defined in the previous chapter 2.4.2). Although this method has been proven to be efficient to determine the melting or enthalpy of fusion of a given material, the resulting specific heat profile as a function of the temperature depends on several parameters of the experiment [2][3][4] being the mass of the sample used for the experiment and of the heating rate of the DSC. To avoid such dependency, Günther et al. [2] used a so-called DSC isothermal step method that consists in realizing a series of measurements for different temperature steps rather than on a temperature range. This method was shown to be particularly efficient for low conductivity materials similar to the PCM used here and is used to characterize four different PCM of interest.

Once the different PCM are characterized, the use-case can be modeled and simulated to assess its performances and optimize the cavity geometry in which the PCM composite is enclosed. For that purpose, a series of parametric analysis is realized to determine the influence of the geometry enclosure, the filler porosity and material as well as the optimum PCM among the four already characterized. The results of the parametric study are compared with the performances obtained with the current solution in terms of both end CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE of mission temperature and mass. To confirm the results obtained during the previous parametric analysis, a series of experiments is run. These experiments aim at reproducing as accurately as possible the conditions of the use-case. To do so, eight samples with different enclosure geometries are manufactured through 3D-printing and filled with Rubitherm RT90HC.

In the first part of this chapter, four different PCM with four different melting temperatures are characterized and their specific heat capacities are determined as a function of temperature. In the second part, the parametric study of the PCM enclosure is realized to determine its optimum geometry. Finally, the experiment is performed to assess the simulation accuracy and validate the system viability as thermal management solution.

Phase Change Material Characterization

In order to assess the viability of PCM based system for the use-case described in introduction, it is first necessary to characterize one or more PCM that may be adapted to it. For that purpose a series of measurements are realized on four different PCM which characteristic are summarized in table 3.1. Their melting temperatures are chosen to be in range determined for the use-case and having the highest latent heat as possible. The selected PCM are tested using two methods: DSC dynamic and isothermal step methods to avoid influence of the measurement parameters in the melting range while saving time in the solid and liquid zone as explained later. Unlike in the previous chapter, it should be noted that all experiments described in this chapter has been performed at IEMN laboratory.

Differential Scanning Calorimetry Dynamic Method

The first series of measurements consist in realizing a standard DSC on the different PCM.

As ramp: one is left empty and serves as a reference while the other is filled with the sample to be characterized. One of the issues linked with such method was pointed out by Günther et al. [2] who studied such method for the measurement of the melting temperature and latent heat of fusion of low thermal conductivity polymers similar to the PCM studied here.

The problem of the DSC method lies in the way the heat is measured: it is done thanks to a thermocouple placed underneath the crucible made of aluminum. When the chamber in which both crucibles are placed is heated, the temperature of the crucible rise up quickly following the one of the chamber but if this crucible is filled with a low conductivity material, it can take some time before the temperature reaches a similar level inside it. This difference in thermal conductivity creates a temperature gradient inside the system (crucible + sample) and the measurement is biased. For stronger heating ramp the effect is even stronger as the delay in temperature setting accumulate over time. A similar conclusion can be made for heavier samples which store more heat. This effect is similar to what was observed in part 1.

regarding Thermal Lagging or Local Thermal

Non-Equilibrium (LTNE).

To confirm Gunther's observations, it was decided to test all four PCM with different heating rates H r (i.e. 1 K•min -1 , 5 K•min -1 and 10 K•min -1 ) and sample mass (reported in table 3.1) leading to a total of six measurements per PCM. We depict in figure 3 of the sample can be observed on the heat flow resulting curve corresponding to the two peaks (respectively downward and upward) and depicting heat absorption of heat release.

A further study of figure 3.2 allows to see the influence of the sample mass on the heat flow response: for the heaviest sample, the energy is absorbed on a wider temperature range showing the effect of thermal lagging. In that situation the actual temperature of the PCM is different from the one measured underneath the crucible and the effect increase for higher mass. To derive the specific heat capacity from the heat flow and temperature curves, the following relation is used:

C p (T) = 1 m ∂ ∂T t 0 q (T(t )) d t , (3.1)
with m the sample mass, q the heat flux, T the temperature and t the time.

The resulting profiles of specific heat capacity as a function of temperature are represented in figure 3.3 (similar curves for the three other PCM are found in Appendix A. 1.4) in which the combined effect of heating rate and mass appears clearly. For heavier samples at a given heating rate the melting of the PCM seems to occur on a larger temperature range although this property is a material characteristic. The same conclusion arise from the study of different heating rates for a given sample mass. Although, choosing lower heating rate and mass seems to be a way to draw near the actual specific heat capacity profile, issues can arise. First as seen in figure 3.3, for lower heating rate, measurement uncertainty arise, especially outside the melting range which is due to an increased signalto-noise ratio increasing the error on the estimation of the heat flux [2]. Decreasing the mass of the sample could also be seen as a solution but it should be reminded that most PCM are not pure material which in the case of very light samples could lead to a slight change of composition from one sample to another. These examples shows that another method should be found to characterize the different PCM.

Differential Scanning Calorimetry Isothermal Step Method

In the previous paragraph, it was shown that DSC dynamic method fails to actually depict the specific heat capacity profile due to the influence of the heating ramp and sample mass. As it is nearly impossible to determine the optimum parameter another method needs to be used. Ghünter et al. proposed a different method called the DSC Isothermal

Step method. Whereas a sample is subjected to a temperature ramp while using the DSC dynamic method, temperature profile is different is the DSC Isothermal Step method as it is subjected to a succession of isothermal steps. From one temperature step to another, temperature is slowly increased and then kept constant until the heat flow falls down to near zero. This allows for the temperature to stabilize in the entire system (crucible + sample) avoiding any disturbance by the temperature gradient. The value of the temperature step, called the resolution, can be tuned to refine the measurement and obtain a greater number of points. Although it is tempting to increase the resolution (or decreasing the temperature step) to obtain as many points as possible, it should noted that, first, this is time consuming and secondly, for very high resolution values, the energy absorbed for a given temperature step vanishes which, in general, decrease the precision of the measurement. For each temperature step, a heat flux peak (see figure 3.4a) is observed, which corresponds to the energy absorbed at this temperature step (h i ) (see figure 3.4b) (similar curves for the three other PCM are found in Appendix A.1.4) which is calculated using:

h i = 1 m t i t i -1 q (T(t )) d t , (3.2)
with t i the ending time of the i th temperature step. In this study, in order to save time while providing enough time for the energy to be stored by the system, irregular time step were set. This is seen in figure 3.4a as time step are larger around fusion temperature where the more energy is to be stored.

From the different value of energy absorbed at each temperature, the latent heat of fusion can be estimated by simply summing the value of enthalpy of the highest peaks.

Here we estimated the enthalpy of fusion to be 188 J•g -1 to be compared to the 180 J•g -1 announced by the manufacturer, which tends to validate the method. From those same values, the specific heat capacity can be determined at the different temperature steps using the following relation:

C p (T i ) = h i -h i -1 T step , (3.3)
with T i the temperature of the i th temperature step, h i the energy absorbed by the sample at at the i th temperature step and T step the resolution. The values of specific heat derived from the DSC isothermal step method are reported on figure 3.5 to be compared with the results obtained with the DSC dynamic method. As seen, on this figure, the DSC

Isothermal

Step method is only realized on a relatively small temperature range as it is more time consuming than the DSC dynamic method. Comparing the two methods, the CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE DSC Isothermal step method gives a tighter and higher peak in the melting range than the dynamic method but on a smaller range and reduced point number. While the values of specific heat capacity derived from the dynamic seems erroneous, at least in the melting range, the values derived from the isothermal step method can simply not be used as an apparent specific heat, or at least not as so.

DSC Methods Combination and Effective Specific Heat Capacity

Modeling

In the two previous paragraphs, it was shown that while the DSC dynamic method was too dependent on the experiment parameters to be accurate inside the melting range, it gives repeatable and accurate specific heat values outside this range when the material is either completely solid or liquid. On the other hand, the DSC isothermal step method precision only depends on the value of the temperature step, that should be low enough to maximise the number of data point but high enough for the measurement precision to be suitable. With a temperature step of 0.5 • C, we showed that a good approximate of the latent heat of fusion and specific heat could be obtained inside the melting range.

Outside this range, only few data points are available due to increased measurement time compared to dynamic method.

To derive an accurate representation the apparent specific heat from the two different sets, it is interesting to combine their advantages: a quick and reliable specific heat measurement outside the melting range using the dynamic method, and a slower, discrete but more precise measurement of it inside the melting range using the isothermal step method. Combining those two requires to find a suitable general expression of the specific heat that could be fitted to the experimental results. Only a few studies can be found in the literature regarding such function. One of them was defined by Yang et al. [5],

who represented the enthalpy of the material a a function of temperature using a Sigmoid function to smooth it around the melting temperature, avoiding a brutal step jump that cause oscillations because of the difficulty in evaluating the capacitance matrix. Based on the experimental results from the different DSC methods, another expression was thought to be adapted which is derived from the Pseudo-Voigt function, a combination of a Gaussian and a Lorentzian functions which general expression is given as so: The exact same method is used on the three Rubitherm PCM reported in table 3.1.

C p (T) =              a s T + b s + φ π 1 (T-T f ) 2 + α 2 , T ≤ T f , a l T + b l + ψ γ 2 exp - T-T f γ 2 2 , T > T f , (3.
Their resulting specific heat capacity profiles are also depicted in figure 3.5b and compared to the one obtained for the PureTemp 108. The different melting temperatures are noticeable and confirm the manufacturer values. As it was done for PureTemp 108, we notice similar tight melting peak denoting a high PCM quality.

Use-case Parametric Study

In the previous chapter and paragraph, different modeling tools (i.e. ETC model for architected enhancers and apparent specific heat model for four different PCM) were developed. The next paragraph intends to use those tools to assess the possibility of using such enhanced PCM on a use-case of industrial interest through a series of simulations. In addition, the study should provide ways to optimize the PCM based system for the given application.

Use-case Definition and Physical Model

To better illustrate the interest of the above mentioned models, we use it to simulate the behavior of a high power electronic thermal management device. This device aims at limiting the temperature rise of four dies developing each 100 W during a 30 s time span in relatively tough environment at a temperature of 80 • C. The PCM based system is made of a composite composed of a PCM infused aluminum lattice structure which is enclosed in an aluminum package (AS7 alloy). Based on the results of the previous chapter (see part 2.2.3), a FCC lattice is selected with a porosity ranging from 0.7 to 0.9, which allows sufficiently high thermal conductivity (at least 15 times the one of the pure PCM) while leaving a large part of the volume for the PCM. Regarding the overall shape of the system, it is cylindrical with a fixed radius R of 60 mm and total height H varying from 9 to 11 mm (see figure 3.6a). The PCM based composite itself is enclosed in a cavity inside the system that is described in terms of the thickness of the different walls surrounding it (see

figure 3.6b).
The optimization itself should determine the optimum value of different parameters that minimize the end-of-mission die temperature. We list below the different parameters CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE to be optimized:

• System overall thickness H,

• Filler porosity ε,

• Top cover thickness h t ,

• Bottom cover thickness h b ,

• Side wall thickness h s ,

• PCM choice between the four studied in part 3.2.

To determine the end-of-mission temperature of the dies (i.e. after the 30 s time span)

we need to solve the heat equation on the above mentioned device using COMSOL MULTI-PHYSICS™. For that, we divide the overall system called Ω into two different sub-domains Ω 1 and Ω 2 respectively representing the enclosure made of bare aluminum (AS7 alloy) (see grey part in figure 3.6b) and the enhanced PCM (see blue part on figure 3.6b). In Ω 1 the heat equations in a solid is solved leading to:

k s ∇ 2 T = ρ s C p,s ∂T ∂t , (3.5) with:
• T: Temperature [K],

• t : Time [s],

•

k s = 137 W•K -1 •m -1 : Aluminum thermal conductivity, • ρ s = 2542 kg•m -3 : Aluminum density, • C p,s = 900 J•K -1 •kg -1 : Aluminum specific heat capacity.
In Ω 2 , it is necessary to account for both fusion of the PCM and for the presence of both the PCM and the lattice structure. Dealing with the former, we used a fixed-grid method (or so-called enthalpy method) defined by Voller [6] and specific heat profile extracted from DSC measurement in the previous paragraph. To model the behavior of the PCM infused porous network we use a simple One-Temperature model that assumes Local Thermal Equilibrium (LTE) between the PCM and the porous network which is supposed to be precise enough in this case but should be confirmed experimentally in the next paragraph. In that case the heat equation can be written:

K eff ∇ 2 T = (C p, f ρ f + C p,s ρ s (1 -) ∂T ∂t . ( 3.6)

CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE

To account for both heat sources and influence of the environment, the following boundary conditions are set. On Σ 1 , which corresponds to the surface in contact with the dies, the heat flux is fixed as follows:

-

k s ∇T.n = P 0 S 0 , (3.7) with:
• P 0 =100 W: Power delivered by a single die,

• S 0 : Area of a die in contact with the PCM based device.

The surface Σ 2 (in blue on figure 3.6a), which corresponds to the remaining outer surface of the device is subjected to a general convective boundary condition yielding:

-

k s ∇T.n = h(T -T ext ), (3.8)
with:

• h = 20 W•K -1 •m -2 : Convective exchange coefficient, • T ext = 80 • C: External environment temperature.
To better illustrate the performance of such device, it will be compared to the current solution used on such issue which consists in a bare aluminum (Al 6061 alloy) disk, similar to the PCM based system. We report in table 3.2 the thermal properties of the different materials used in the simulations.

Material Thermal Conductivity [W•m -1 •K -1 ] Density [kg•m -1 ] Specific Heat [J•K -1 •kg -1 ]
Al Table 3.2: Thermal properties of the different materials used in simulation.

Geometry Parametric Analysis Height, porosity and top cover thickness influence

We start the parametric analysis of the device by testing the influence of device height, top cover thickness and filler porosity while both bottom cover and side walls thicknesses are fixed at 1 mm. Three different device heights were selected (9, 10 and 11 mm) as it corresponds to a reasonable allocated space for such application. The top cover thickness is chosen so that its maximum value for a given device height leaves a 2.5 mm height space for the enhance PCM while its minimum value will be 1 mm. This value was chosen as it corresponds to the smallest pore size that can be manufactured for the level of porosity required in this application. Lastly, filler porosity ranges between 0.7 and 0.9 which is a good compromise between high ETC and high fraction of PCM in the composite.

We depict on figure 3.7a the mean temperature on the die surface in contact with the heat spreader at the end of the mission for all the above mentioned configurations as well as the one obtained for bare aluminum devices of different heights (denoted "Bulk"). The latter serve as reference as it is an usual heat spreader used in many industries.

We can first see the influence of the height of the device: a higher device provide an increase mass that helps store additional heat allowing lower temperature at the end of the mission. Though it is promising, the additional mass is, in most industries, to be avoided. This can be seen for Bulk configurations as well as for PCM devices no matter the configuration.

If we now look at the influence of the top cover thickness we see that no matter the device height or filler porosity an optimum of this value exists. We especially notice that the configurations containing the most PCM are not the ones that performed the best. This could be explained by analyzing the melting ratio ( f l ) at the end of the mission for the different configurations depicted in figure 3.7b. The melting ratio is defined as the frac-CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE tion of PCM melted in a given sample at a given time and evaluated using equation 1.27 at each point of the discrete Finite Element Analysis (FEA) network and then integrated over the entire structure. We notice that for low values of the top cover thickness, up to 35% of the total amount of PCM is still un-melted at the end of the mission. This means that some latent heat that would have been able to absorb heat is wasted. This is due to the fact that heat cannot flow easily through the relatively bad conductor that the enhanced PCM is. Increasing the top cover thickness seems to enhance the ability of the device to efficiently spread heat away from the die, heat that can be stored in the form of latent heat once it reaches and melts PCM. On the other hand, for the highest values of top cover thickness, all the PCM has melted before the end of the mission. This is detrimental in terms of thermal management as once melting has occurred, enhanced PCM remain a relatively bad thermal conductor avoiding an efficient heat spreading from the dies which increases their temperature. The optimum value of top cover thickness is actually the optimum between an efficient heat spreading to avoid heat concentration around the die while providing a decent amount of PCM to absorb as much heat as possible.

The influence of the porosity is harder to determine. Comparing the optimum value of temperature for a given device height we notice a slight change from one porosity to another while a slight shift in terms of optimum top cover thickness is observed. The change in optimum top cover thickness depending on the porosity is understandable as for lower porosity the enhanced PCM has a higher ETC allowing better heat flow which necessitates thinner top cover for similar results. This is well illustrated if one looks at the extreme values of top cover thicknesses: for very thin top cover (i.e 1 mm), the higher ETC of the composite with a 0.7 porosity allows for better performance than composites with higher porosity. The reverse is also true when looking at high value of top cover: in this case higher porosity allows for more space to store PCM increasing the device performances.

Looking at the optimum geometries and comparing them to the temperature obtained with bare aluminum devices for the different values of overall thickness, we observe a gain of respectively:

• 10.4 • C and 84.7 g (about 30% gain) for a 9 mm thick device, 4 mm thick top cover and a composite porosity of 0.9,

• 10.4 • C and 95.0 g (about 31% gain) for a 10 mm thick device, 4.5 mm thick top cover and a composite porosity of 0.9,

• 10.3 • C and 80.6 g (about 24% gain) for a 11 mm thick device, 4.5 mm thick top cover and a composite porosity of 0.8.

Bottom cover and side walls thickness influence

Similar simulations are run to understand the influence of the side walls and bottom cover thicknesses. To do so, we fix porosity at 0.9, bottom cover thickness at 1 mm and device height at 9 mm and vary both top cover thickness from 1 to 5.5 mm and side walls thickness from 1 to 4 mm. Results of those simulations are depicted n figure 3.8. In this figure, we see that side wall thickness only have small influence on the performance of the device which could be explained by the fact that those walls are far from the die compare to the top cover and they participate less in heat spreading. Increasing their thickness only reduce the space available for the composite, hence reducing the amount of heat that can be stored through latent heat.

Another set of simulation is ran this time only varying bottom and top cover thicknesses. Four cases are studied for four different combination of top and cover thickness which are respectively set at: to the top cover, the bottom one is positioned far from the dies and does not participate, at least at the beginning of the mission, to the heat spreading. So the bottom cover is to be kept as thin as possible.

PCM Choice

We finally compare the behavior of the four PCM previously mentioned to see how the melting temperature and latent heat of fusion could affect the dies temperature profile and finally their temperature at the end of the mission. For that, we consider a device similar to the one in figure 3.6a with a height of 10 mm, a top cover thickness of 4.5 mm, side wall and bottom cover thicknesses are fixed at 1 mm and porosity at 0.8. We depict the variation of the dies temperature profiles in figure 3.9. In the figure, we can actually observe the melting of most PCM occurring as the temperature profile's slope is modified and reduced during the melting phase before increasing again when melting is over in the entire enclosure. If we now compare the different PCM we notice that Rubitherm RT 80HC is the one that performs the worst despite its higher latent heat of fusion. This actually makes sens as its melting temperature is the one of the simulation environment hence no melting happens during the mission in this particular case. Comparing the other PCM we notice that Rubitherm RT 90HC has a slightly lower latent heat of fusion than Rubitherm RT 100HC which explains why the latter outperforms the former. Its latent heat of fusion is also lower than PureTemp 108 but as its melting temperature is lower, its melting starts earlier which seems to be favorable, at least for relatively short and brutal heating as the one studied here.

C)

Lattice : ǫ = 0.7, H = 9mm Lattice : ǫ = 0.7, H = 10mm Lattice : ǫ = 0.7, H = 11mm Lattice : ǫ = 0.8, H = 9mm Lattice : ǫ = 0.8, H = 10mm Lattice : ǫ = 0.8, H = 11mm Lattice : ǫ = 0.9, H = 9mm Lattice : ǫ = 0.9, H = 10mm Lattice : ǫ = 0.9, H = 11mm

(a)

1 2 3 4 5 6 7
Top cover thickness (mm) A first observation can be made regarding a similar optimum value of top cover thickness for any case of overall thickness and porosity but the profiles are different. Then, one can notice the stronger influence of the porosity for the lowest value of top cover thickness that could be related to the value of the melting temperature: based on Fourier law, any decrease in ETC automatically increases the temperature gradient in the system and consequently in the composite. If the thermal gradient is higher, it usually means that more

PCM will be at a temperature below the melting temperature at which most heat is stored leading to a higher temperature at the heat source. Regarding the optimum value of top cover thickness, a shift toward higher values is noticed compared to similar simulations with Rubitherm RT 90HC. The explanation can be found in figure 3.10b, in which are displayed the melting ratio at the end-of-mission for the different configurations. Contrary to the same simulation with Rubitherm RT 90HC, no configuration has seen a complete melt of the PCM which means that the optimum occurs before it. The fact that Rubitherm CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE RT 100HC performs so well despite the melting process not being over might actually be related to a its higher latent heat of fusion compared to Rubitherm RT 90HC.

Comparison with different types of filler

Discussion

In the previous paragraphs, we were able to illustrate the use of an ETC model for lattice structure used as PCM filler. We show how the above mentioned lattice structures could be used to manufacture a thermal management device for high power electronics in tough environments. We demonstrate that for the device to actually outperform current solutions like bare aluminum heat spreader, a refined composite enclosure optimization was required. In that sens, we studied the influence of six parameters:

• Device height,

• Top cover thickness,

• Bottom cover thickness,

• Side walls thickness,

• Filler porosity,

• PCM choice.

The study showed that increasing the device height would reduce the electronics temperature but bring additional mass which in certain cases could be troublesome. Then, it was demonstrated that an optimum of the top cover thickness exists which allows sufficient heat spreading while allocating enough space to the PCM. On the other hand both bottom cover and side wall thicknesses should be kept to minimum as they do not participate to the heat spreading but could reduce the amount of PCM stored in the device. It was shown that the filler porosity (within the studied range) was not a key parameter to be optimized as performance were equivalent for all porosity values tested here.

From a quantitative point of view, the parametric optimization allowed an average 10 • C and a 30% mass gain compared to the current solution made of a bare aluminium disk. The temperature gain is to be compared to the total thermal heating which represents about 56 and as shown by Yang et al. [7], its influence highly depends on the orientation of the heat source. In this use-case, the heat sources are on top of the sample which should limit the influence of convection to minimum, making the assumption of absence of convection correct. On the other hand, even if LTE assumption has been proven correct in a variety of cases [8][9][10][11][12] it might not be the case here depending on a parameter that has not been discussed in the previous paragraph. Actually, for a better integration and to avoid any contact resistance between the package and the composite itself, the final device will be manufactured through 3D-printing. One issue related with additive manufacturing is the limitation in terms of size reduction which is due to the metallic powder (grain size, material...) and laser used for the powder sintering (see part 2.4.1). To maintain the desired level of porosity, one must maintain a sufficiently high pore size leading to a lower specific area which has been proven to favor LTNE (see part 1.3.2). Assessing the validity of the LTE assumption could be done through Direct Numerical Simulation (DNS), although it would require a huge computing power. Another way to do so is through experiments which is the subject of the next part.

Experimental Validation of the Model Implementation

Samples description

Sample design

To validate the different assumptions made in the previous paragraph as well as the overall simulations assessing the PCM based devices, it was decided to experimentally test several such devices. Similarly to the previous samples (see part 2.4.1) such devices are produced through additive manufacturing, which is the only way for the intricate topology of the lattice to be perfectly integrated into the aluminum package. To avoid any possible additional thermal resistance between the lattice and the package, that can be consequence of bad thermal contact, the entire package is manufactured as one block which is seen on the tomography images (figure 3.11) taken on Sample #3.

The different samples geometrical parameters are defined in table 3.3. Those samples were designed to experimentally assess the influence of the following parameters on the end-of-mission die temperature profile:

Sample printing & filling

Once designed, the samples are produced through additive manufacturing, using the Selective Laser Melting (SLM) method that consists in building a given structure from a Computer Aided Design (CAD) file by locally melting parts of powder layers thanks to a laser beam [13]. The building process starts with a powder layer on a substrate that is melted by a high energy density laser beam at selected location defined in the CAD file. Once the first layer has been melted at the right locations another layer is added on top of it and melted at selected locations again. The process is repeated until the entire structure is printed. The remaining powder left un-melted, needs to be removed from the inside of the package. For that reason, eight holes (see figure 3.11) had to be added on the edges of the top cover. As seen in figure 3.11b, powder extraction was satisfactory as only few powder grains remain.

All samples were filled with Rubitherm RT 90HC which filling was performed at 125 • C; although this is above the maximum operation temperature recommended by the manufacturer, it ensures low viscosity during the filling and optimum filling ratio. The filing is performed through the same holes used to extract powder which are then sealed using sealing plug made of steel ball inserted in a steel tube which inner diameter is slightly smaller than the one of the ball. Once inserted the plug is maintained by the pressure applied by ball on the cylinder and holes wall. These sealing plugs requires a sufficiently thick and robust baseplate to be sealed in. For that reason, and in respect to the sealing plugs dimensions, a cube is added underneath each hole of the sample (see figure 3.11b).

They locally reduced the PCM allocated space but allows for robust sealing and avoid any leaking.

In addition to the geometry related parameters, the roughness of the sample surface is also a key point that may influence experimental results. Actually, to heat the samples up a heating resistance (described in the next paragraph) with a contact surface area similar to the one used in the simulation is used. In using such heating method, the resulting heat flux delivered by the resistance to the device greatly depends on the surface state.

Depending on the sample, a noticeable variability of the surface state is noticed which could be seen in the experimental results later.

Experimental Setup Description

The experimental setup needs to reproduce the conditions in which the use-case should perform while ensuring precise sample surface temperature measurement and reliable thermal contact between the heat elements and the sample.

To do so, four heating resistances with a maximum delivered power of 140 W are used (reference : NIKKHOM -TO247 140W High Power Resistor). To ensure good thermal contact they are not directly put in contact with the sample but a layer of thermal grease with a thermal conductivity of 4 W•K -1 •m -1 is added. The layer thickness should be sufficient

Experimental Setup Calibration

CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE

The second step of the calibration is to account for heat loss. Actually, while the entirety of the energy enters the sample in the simulation of part 3.3.2, some may not in the experiment. Looking back at figure 3.13, one can notice the large copper wire section which is about 2.5 mm 2 used to allow for safe powering as current could be as high as 25 A.

Although it is necessary, the thermally conductive copper drains a part of the power delivered to the dies. To model it, the total power delivered to the heating element is reduced by 15%. In addition, natural convection is considered on top of the die stack-up, with a convective exchange coefficient of 400 W•m -2 •K -1 . It should be noted that this coefficient does not represent actual convection but is rather a way to model heat loss by convection and conduction. These values were found to fit experimental data as shown in figure 3.15 but were not actually determined experimentally. Figure 3.14 summarizes the different assumptions as well as the final geometry of the model used in the simulation.

On the experimental side, the temperature is measured at four different locations as seen in figure 3.13. Although the four thermocouples are supposed to give similar profiles, they tend to be slightly different because of the uncertainty on the actual of both thermocouple and die locations. Actually, because of the high power density the temperature gradient, especially close to the die, is relatively high and the resulting temperature profile could be different for two relatively close locations. To limit the uncertainty, an average of the temperatures measured with the four thermocouples is used to be compared to the simulation results.

In figure 3.15 are reported both simulation and experimental data for the three different bare aluminum samples. Using the assumptions described above, a maximum 5% deviation is obtained between simulation and experiment which validates both simulation's boundary conditions and experimental setup. Due to the reduced power that actually enters the dies, the end-of-mission temperature is lower than the one predicted in the previous paragraph and depicted in figure 3.7 which means that results from these experiments can only be used to verify the validity of the assumptions and, if necessary, correcting them.

Experimental Results and Discussion

Now that the simulation is calibrated, both experiment and simulation can be run on the PCM based devices to assess the validity of the different assumptions previously made:

LTE and absence of fusion-induced convection.

Influence of the pore size

To assess the necessity to consider LTNE or not, the first three samples to be studied are #3, #5 and #8 which only differ by their pore size with fixed porosity, top cover thickness and overall thickness. As reported in table 3.4, and due to constriction, the pore size has an influence on the ETC and variation of about 15% are observed in the vertical direction between sample #3 and #8. Using these properties with a LTE assumptions and the same set of equations than in part 3.3.1 gives the results shown in figure 3.16a. On that graph, we notice that despite the difference in ETC, the results are supposed to be similar. These could be compared to the experimental results depicted on the same graph and display high deviation (up to 30%), with the larger deviation obtained for the biggest pore size.

To understand the deviation, we observe the different experimental temperature profiles which shows that the bigger the pore size, the higher the deviation. Apart from the ETC, the only parameter varying between those devices is the pore size hence the specific area. Data in table 3.4 show that specific area in sample #3 is almost three times as high as in #8 which can be directly related to LTNE through equation 1. 48. In order to fit experimental data, simulations are run again but this time with a LTNE assumption which requires equation 3.6 to be replaced by the following set of equations:

CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE              R x y 0 0 0 R x y 0 0 0 R z k s (1 -ε)∇ 2 T s = (1 -ε)ρ s C p,s ∂T s ∂t + h(T s -T f ), k f ε∇ 2 T f = ερ f C p, f (T f ) ∂T s ∂t -h(T s -T f ), (3.10)
with: R x y = G s,x y and R z = G s,z 1 +

2R c P s G s,z (1-ε)k s N -1
and h an integral heat transfer coefficient for the heat conduction at the fluid-solid interface that is to be determined for each

sample. An estimation of such parameter (see equation 1.49) was done by Quintard and

Whitaker [14] on an array of cylinder by solving a so-called closure problem. Although the FCC cell displays a different topology, such expression was used as a baseline to approach the actual value of h that was then modify to better fit experimental data. We did so for sample #3, #5 and #8 with respective value of 9.05

× 10 5 W•m -3 •K -1 , 5.6 × 10 5 W•m -3 •K -1
and 2 × 10 5 W•m -3 •K -1 which is coherent as h increases as the pore size decreases almost following 1/P s a proportionality law. Modifying the LTE assumptions showed great results as deviation is maintained under 5% for the entirety of the mission with an increased difference at the end of the mission that might be due to a lower filling ratio than expected.

Actually, due to the PCM expansion at the phase transition, the PCM enclosure cannot be full to avoid over-pressure at the liquid state which could be detrimental to the mechanical integrity of the device. Using those assumptions, we notice despite all other parameters being fixed, sample #3 with the lowest pore size outperforms the two others due to its higher specific area that tends to limit thermal lagging.

Influence of the porosity

We just experimentally proved that pore size has actually an influence on the Latent Heat Thermal Energy Storage System (LHTESS) overall performance due to thermal lagging which tends to worsen the devices performances and could be accounted for considering a LTNE assumption. We continue the study, this time considering the influence of porosity only on the performances of the LHTESS. To do so, samples #4, #5 and #6, with respective porosity of 0.7, 0.8 and 0.9, are tested and compared. Similarly to what was done in the previous paragraph, the three samples are simulated using both LTE and LTNE assumptions, the latter requiring the determination of a coefficient h, to see if one of them could fit experimental data.

Starting with the LTE assumption (see figure 3.17 uncertainty on both thermocouple and die position that might affect the measurements, it is considered satisfying.

Although we overestimated the performance of the LHTESS in the previous paragraph, it is interesting to notice that it was already predicted that porosity would have minimal influence on the end-of-mission temperature for near optimum top cover thickness. The combined effect of porosity and thermal lagging seems to have made this difference even smaller as samples #3 and #4 temperature profile are almost identical.

Influence of the overall package thickness

Comparing samples #1, #2 and #5 (see figure 3.18), we want to assess the influence of the overall package height on the temperature profile of the die. Similarly to the two previous paragraphs, we start the study by comparing the results derived from simulations using the LTE assumption and the experimental data. Also the simulations are able to display the influence of the overall height, hence mass, of the system they once again fail to fit experimental data as an up to 18% deviation is observed.

To enhance the simulations, a LTNE assumption is used with respective h values of 6.32 × 10 5 , 3.58 × 10 5 and 5.6 × 10 5 W•m -3 •K -1 reaching satisfying level of deviation (<5%).

The relatively low value of h for both samples #1 and #2 cannot be explained by Quintard and Whitaker model but should be related to the coefficient R defined in part 1.3.2 (see equation 1.45). This coefficient, characteristic of the local thermal equilibrium, is to be minimized if LTE is to be considered and is inversely proportional to the system size, hence the device thickness. In addition, it should be noted that the delayed temperature increased observable for sample #1 may be related to its surface state. Actually, due to the way it is manufactured, the surface state might be relatively rough. This creates an additional thermal resistance at the interface between the heating element and the sample itself, increasing the heating element temperature at the beginning of the experiment and reducing the sample surface temperature. This explains the unexpected low end-ofmission temperature compared to sample #2 and #5 as well as the higher relative error.

Influence of the top cover thickness

Finally, the influence of the top cover thickness is assessed comparing samples #5 and #7.

Those samples differ in terms of top cover thickness and in terms of pore size. Based on previous observations and on the results of the previous simulations both using LTE and LTNE a high difference is expected.

Once again, even though LTE simulations correctly predict a higher end-of-mission temperature for sample #7, they fail to fit experimental data, reaching an up 25% deviation (see figure 3.19b) and LTNE is necessary. Using such assumption and respective h values of 5.6 × 10 5 and 2.37 × 10 5 W•m -3 •K -1 we are able to fit experimental measurement with a maximum 5% error. It should be noticed that once again discrepancies are observed at the beginning of the measurement for sample #7 due to its rougher surface creating additional thermal resistance.

Overall considerations

Before going further into the performance assessment of such solutions, it is important to mention that while the use of LTNE helped fit experimental data, the integral heat transfer coefficient (h) determination is prone to errors. Actually, the different coefficients determined here accounts for Thermal Lagging but they also take into account experimental errors (thermal resistances and thermocouples misplacements, surface roughness leading to additional thermal resistances...). For that reason, we thought it necessary to com-Sample h (W•m -3 pare those values with those predicted by Quintard and Whitaker [14] for a cylindrical system using the following equation:

h = 8k f 2 P 2 s (1 -) k f (1-)k s + 1 -3 -4 ln . ( 3

Lattice/PCM composite based device performance assessment

Experimental study of the eight different samples showed that LTNE assumptions has to be considered, for the simulation to accurately predict the temperature evolution on the sample surface. Adding thermal lagging has been proven to weaken the device performances and the enhancements found in part 3.3.2 may be impacted. To assess the impact of thermal lagging, we perform again the simulation done in part 3. As a consequence, it could be concluded that aluminum lattice based filler itself might not solve the issue and further topological enhancement is necessary. As mentioned in chapter 1, Cascaded storage, which consists in using different types of filler or graded filler in a single system to enhance heat conduction and storage, could be used to further enhance the system storage efficiency which will be studied in the next chapter.

Conclusion

This chapter aimed at demonstrating how the modelling techniques derived from the literature and developed in the previous chapter could be implemented in a simulation software (COMSOL MULTIPHYSICS™) in order to assess and optimize the performances of a PCM based device for the thermal management of power electronics.

One way to model PCM fusion is the use of an effective specific heat capacity that would vary with the temperature, especially close to the melting temperature [15]. In the first part of the chapter four PCM are characterized in that sense: their specific heat capacity is measured through two types of DSC method. The first one, called the dynamic DSC subjects the sample to be characterized to a constant temperature rate and measures the energy absorbed. The issue with such method is the bias brought by the difference of thermal conductivity between the crucible in which the sample is placed and the sample itself. For heavier samples or faster temperature ramps, this tends to create a temperature gradient in the system that skews the measurement, especially around the melting temperature. To avoid such bias another method had been developed: the isothermal step DSC method that consists in subjecting the sample to a series of isothermal steps at which the absorbed energy is measured. Although this method avoids the measurement bias, it gives a discrete representation of the specific heat capacity on a reduced temperature range. So, it was decided to combine both methods: using the dynamic method far from the melting range and the isothermal step method close to it. Fitting the experimental data to a pseudo-Voigt model, we were able to determine the specific heat profiles for the four PCM.

In the second part of the chapter, the different modeling tools developed upstream are used to assess the performance of a PCM device on a specific use-case. This use-case aims at maintaining four dies each emitting 100 W for 30 s in an environment at 80 • C under 120 • C. For that the space available is cylindrical-shaped with a diameter of 120 mm with an up to 11 mm height. The device itself is made of an aluminum package in which is CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE enclosed an aluminum lattice/PCM composite. Through a parametric study, the geometry of the enclosure and the composite are optimized and it is determined that:

• the top cover of the package, serving as an interface between the composite and the heat source, possesses an optimal thickness allowing an efficient heat spreading through the conductive aluminum packaging while providing enough space for the PCM for the heat to be absorbed,

• all other covers (side walls and bottom cover) should have minimal thickness as they barely participate to the heat spreading, being far from the heat source but reduce the space allocated to the PCM,

• moderate variation of the porosity (from 0.8 to 0.9) seems to have almost no influence on the performance of such system.

To confirm the trends determined in the previous part and validate the different boundary conditions and assumptions, a series of samples were designed, additively manufactured and experimentally tested. Experimental results showed that the simulations made in the previous part were over-simplified which tended to over-estimate the devices performances. Actually, thermal lagging was neglected in the parametric study yet it usually tends to slow down the fusion process which cause slower energy absorption and increased temperature. Using a LTNE assumption instead of the LTE assumption made in the previous part allowed to fit experimental data although it confirmed that device performance were over-estimated. This shows that the enhanced PCM used in this part is not sufficient for this given use-case and needs further improvement. This will be the subject of the next chapter that intends to use cascaded storage as a way to further enhance the system used here.

Introduction

In device design in which is introduced a heat spreader. This additional piece aims at better spread heat away from the die in order for it to be stored in a more efficient way. Adding such spreader obviously reduces the space available for the PCM and also increases the device mass and filler porosity needs to be adapted to avoid both. Some simple spreader design are developed and their performances assessed in order to prove the interest of such hybrid filler in reducing the end-of-mission die temperature.

Although they show decent performance enhancement, the simple design of the spreader may not be optimized and improvement might be made. To do so, a biomimicry technique based on leaf vein modeling is used [4][5][6]. Actually, leaf veins are used, in nature, to diffuse sap from the base of the leaf to its extremities in order to supply the entire leaf. Heat and sap diffusion are similar phenomenon and this technique has already been proven to enhance PCM based devices [6]. In the second part of the chapter, we describe how such leaf vein based structures are generated, how they can be tuned and finally optimized using a genetic algorithm. In order to show the interest of such process, it is first used on a simplified 2D study case and the derived structures are compared to a simple lattice filler.

In the third part of the chapter, the previously developed generative design algorithm is applied to the optimization of a finned spreader for the use-case defined in Chapter 3.

We first describe how the spreader generation defined in the preceding part is adapted to fit the geometry of the use-case system enclosure. Then, for time saving reasons, the generative design algorithm is first adapted to a 2D version of the previous use-case in order to select some pseudo-optimized structures that could be adapted in 3D. After discussing CHAPTER 4. BEYOND HOMOGENEOUS FILLER the changes occurring during the transition from 2D to 3D, the resulting structures are inserted inside the LHTESS to assess the improvements brought by the hybrid filler.

The final part of the chapter discuss some of the perspectives that could further improve the performances of such system as well as the modeling tools to more accurately predict its behavior.

Preliminary study

In where q is the heat flux entering the system, ∆t is the melting time and ∆T the temperature difference between the top and the bottom surface of the system. According to Feng et al., hybrid filler possess a heat transfer coefficient seven times higher than the foam and about 24% higher than the array of fins. Although these results are promising, they need to be nuanced as metal fraction in the system change from 3% to 9% and 11.8% going from foam to array of fins and hybrid filler. Yet the calculation of the heat transfer coefficient is actually related to the time necessary for the PCM to melt hence reducing the amount of PCM and increasing the amount of conductive material obviously will reduce the melting time. However, even if this result could be seen as biased, a seven times increase is significant enough for the solution to be considered viable. the enclosure hence the spreader can be defined by its base area A b that can be tuned by playing on the value of the different parameters. In the remaining of the enclosure the aluminum lattice enhanced PCM is still found but to avoid the bias mentioned above, the porosity of the lattice is adapted to keep the same amount of both aluminum and PCM in the enclosure. To calculate the porosity of the filler, we first consider one of the LHTESS from chapter 3 with a filler of porosity ε 0 and if the total base surface of the enclosure is denoted A 0 = π(Rh w ) 2 then the porosity of the hybrid filler ε h (ε 0 ≤ ε h ≤ 1)is expressed as: spond respectively to the maximum temperature found on the die surface throughout the simulation and the average temperature over the surface. In these conditions, the Fin/No filler configuration is outperformed by the two others due to its inability to efficiently conduct heat from the the tip of the spreader to the bare PCM. The comparison of the lattice and hybrid fillers needs to be more subtle. Comparing the average die temperature, as it was done in chapter 3 is not as relevant here. Looking at the average temperature profiles, the role of the heat spreader seems limited. On the other hand, a 2 • C gain at the end of the mission (i.e. 35 s) is observed for the maximum temperature. It appears that adding a spreader provide a way, even if limited, to avoid hot spots that could be detrimental to the electronics reliability. Although interesting, those results were obtained using a LTE assumption, which has been proven in part 3.4.4 to be inadequate to describe the behavior of such device in these conditions. As a consequence, the simulation is performed again for both lattice and hybrid fillers using a LTNE assumption which results are depicted in figure 4.3b. For that purpose the integral heat transfer coefficient defined in equation 1.49 is considered to be equal to 9.05 × 10 5 W•m -3 •K -1 as determined for the best case in part 3.4.4. Although this value has been proven to be correct for the lattice filler, it might be overestimated, due to the higher porosity in the case of the hybrid filler even if only small influence of the porosity on such parameter was found experimentally (see part 3.4.4). As, expected, we notice, in figure 4.3b, an increased temperature compared with LTE assumptions: as described before thermal lagging slows down the melting process and heat storage causing an increased end-of-mission die temperature. However, CHAPTER 4. BEYOND HOMOGENEOUS FILLER in addition to the better heat spreading brought by the hybrid filler, it seems that it also limits the consequences of thermal lagging as seen in figure 4.3b. Due to the higher average thermal conductivity around the die, thermal gradient is smaller and increases at a lower rate. Looking back at equation 1.44 which mathematically describes thermal lagging, this directly affects the overall behavior of the composite by diminishing one of the terms linked to thermal lagging.

ε h =      ε 0 1- A b A 0 , if A b ≤ (1 -ε 0 )A 0 , 1, otherwise.
We are seeing here, that the optimization of the entire system is not only a question of conductive and heat storing materials proportion but also a question of where each of the material should be placed. Although those results are interesting, we notice that progress can be made by determining the best topology for the spreader. One way to do it would be to realize a parametric study on the fin based heat spreader but that may not be sufficient.

Actually, fin shape possibilities are almost endless and classical parametric optimization techniques does not work well, in general, unless a close-to-optimal structure is already known [7]. To go further, topology optimization method is thought to be a good starting point. Topology optimization was originally theorised by Robert Le Ricolais towards a goal that he defined as: "Portée infinie, poids nul" (i.e. zero weight, infinite span). He was the first one to notice the power of structure which consists in finding "where to put holes" in order to produce light and tough structures. Le Ricolais produced a number of pioneer models aiming at studying structure and his work was further improved by Bendsøe [8; 9] who described the mathematical programming of the method based on Finite Element Analysis (FEA). Although both Le Ricolais and Bendsøe used topology optimization to enhance the mechanical aspect of structures, it has more recently been used for thermal management purposes [4; 6; 7; 10-12]. Different techniques are used in order to achieve topology optimization (SIMP, Generative Design Algorithm (GDA)...) and it was decided to use GDA due to its relative easy implementation and already proved efficiency [4]. In the next part of the chapter, we describe the GDA principle and show how it can be used to solve the issue related to our use-case. can be summarized by the following expression:

2D Generative design algorithm

max s F(u(s), s, t f ), subjected to V(s) = V p , R(s) ≥ R min , (4.3)
which denotes that a so-called objective function F, depending on a vector u representing the degrees of freedom of the problem and a vector s describing the design of the shape to be optimized, is to be maximized at a time t f of the simulation. The optimization of F is to be performed under two constraints: the total volume of conductive material V(s) should be equal to a value set at V p and the minimum thickness R(s) of the branches at any given point should be superior to a value R min for the final structure to be manufactured. The actual form of the vector s depends on the algorithm implemented to realize the optimization.

Although such algorithm can be used to solve 3D problems, we start the study on a simplified 2D problem for time constraints. This preliminary study should help us determine the optimum parameter of the algorithm and see the level of performances that are to be expected. In this preliminary study, we want to use such algorithm to find the optimum heat spreader topology to minimize the temperature of a source, that heats a domain Ω depicted in figure 4 Once the temperature distribution is determined, we can define the objective function using the following expression:

F(u(s), s, t f ) = T fus T(x h , t f ) , (4.5)
with T fus the melting temperature of the PCM used in this case, x h represents all the points of Γ h and T(x h , t f ) the main temperature of Γ h (which length is set at L X /40 in this case) at the time t f of the simulation (which is set at 120 s in this case). If t f is wisely chosen it should ensure that 0 ≤ F ≤ 1. The objective function is used to compare the behavior of heat spreaders like the one in figure 4.5 that have to be generated using a space colonization algorithm. Once a heat spreaders have been generated and their objective function have been evaluated, a genetic algorithm is used to select the best ones and makes them evolve for the objective function to be increased in the next step of the iteration. The combination of those two algorithms is called a generative design algorithm which is described in the next paragraph.

Generative design algorithm principle

Generative Design Algorithm (GDA) is defined as "the process of using algorithms to help explore the variants of a design beyond what is currently possible using the traditional design process" [13]. It is an iterative process that test each of the resulting solution and learn from it. The resulting structures are often far from intuitive. The GDA, studied here, is based on Lohan et al. [10] work and depends on two combined algorithms: a so-called space colonization algorithm that will generate the heat spreader based on a series of parameters that are to be optimized by a genetic algorithm.

Space colonization algorithm

The space colonization algorithm aims at filling a given percentage of domain Ω with a continuous structure aiming at spreading heat in the most efficient way. To do so, and following Lohan et al. [10] or Liu et al. [6], we base our approach on leaf vein mimicking. Actually, leaf veins always grow along the most efficient nutrition/energy transfer path with minimal material consumption and seem particularly adapted to heat transfer. Following this principle, it seems that the growth of leaf vein focuses on two criteria: economy of material and efficient configuration [6].

To perform the leaf vein based network, we follow a method developed by Runions et al. [5]. This is an iterative method using two sets of points: the first one represents the so-called auxin points or target points (in red in figures 4.6 and 4.7) and is denoted S. Those N auxins points represent targets that have to be reached by the network at some point in the iterative process. They are the one directing the network growth. On the other hand are found the source nodes (in green on figures 4.6 and 4.7) that constitute the actual network. At the initial step of the network, a given number N 0 Vein of source nodes are given and constitute the current state of the network denoted v 0 . At each step i of the iterative process, each auxin s is paired with the closest vein node of the network v i forming a subgroup s k ∈ S including all auxin nodes associated with the same vein node v k i . From the node pairing process results N growth (N growth ≤ N auxins ) growth directions corresponding to each subgroup of s k that are defined by: Then the created offsprings will undergo Mutation. This genetic operator randomly modifies up to p Mut × N c chromosomes (N c is the total number of chromosomes). This operator allows for better exploration, as some of the alleles randomly introduced may have never been tried otherwise. Once again an optimum of p Mut is to be found as a too low value would restrict the explored region of the optimization space but a too high value would lead to a loss of heredity as the offsprings would only have few common alleles with the parents.

n k = s k ∈S(v k i ) s k -v k i ||s k -v k i || , k ∈ [1, N i Vein ]. ( 4
The final step of the iteration decides whether to repeat or not the process. Many termination conditions could be found in the literature, also we decided to repeat the process for a given number of generation M: higher values of M usually provides better individuals but also increase the computational time and an optimum value for each problem has to be determined.

The main advantages of the genetic algorithm, which flow chart is represented in figure 4.9, lies in two different points. The first one is its adaptability or flexibility: the mathematical requirements of the genetic algorithm are almost nonexistent as it can be adapted to a wide variety of fitness function as long as the vector s can be represented as relatively small number of parameters of all natures. The second one lies in its exploration to exploitation ratio. In optimization processes, the exploitation refers to the use of the previous solutions to guide the evolution of the next ones while the exploration refers to the size of the region of the optimization that has been explored. Usually, genetic algorithms are able to do both with relatively good performances as compared to heuristic (gradient based) algorithm that are highly dependent on the initial step. As seen in figure 4.10, if the initial solution is close to a local optimum, heuristic algorithm will usually converge to this optimum and miss the global one. Genetic algorithm are able to avoid such issue [18] by using a larger initial population that has a higher chance of reaching a point of the optimization space closer to the global optimum.

Phenotype and spreader generation

To combine both space colonization and genetic algorithm, a suitable expression of vector s that will describe the geometry of the spreader and will be optimized, needs to be found. This vector must have relatively low dimensions, for the genetic algorithm to be efficient while precisely describe source and auxin node distribution as well as branch radius and give as much freedom as possible in terms of geometry to avoid any constraint in the optimization process.

In figures 4.6 and 4.7, the number of auxins and source nodes is relatively low which is not the case in general. While it would be possible to detail the coordinates of each source or auxin nodes in that particular example, such method would rapidly increase the dimension of s which is detrimental to the efficiency of the genetic algorithm. Instead, the structure is described using a series of eight parameters:

• D: spatial step size between each node of the network. A higher value of D usually gives more simple or coarser final design, while a smaller value results in a more refined one.

• r min : radius of the extremity nodes which is the smallest value found in the network.

• r max : radius of the source nodes, that corresponds to the maximum value found in the network.

• N sources : Number of source nodes, horizontally aligned along Γ h . This alignment is arbitrary chosen to promote multiple branch growing at the beginning of the space colonization step.

• C Dist : this parameter pick a single distribution among a list of predefined distributions consisting in concentric half circles oriented in different directions (see fig- • N circle : Number of concentric half circles for the chosen distribution.

• N auxins : Number of auxins per half circle.

• α: a parameter used to tweak the shape of a given half circle to reach region of the domain Ω that could not be reached otherwise. (g) α is modified to extend the colonization space, (h) η x = 0.6 and (i) η y = 0.6.

• η x , η y : scaling ratio in x and y direction to reduce the space allocated for space colonization.

To better understand the influence of each parameter on the resulting heat spreader, we use an arbitrary reference spreader with the following parameters that will be modified later on: L x = 20 , D = L x /120, r min = L x /160, r max = L x /40, N sources = 7, C Dist is set for the concentric disk center to be the center of Γ h , N circle = 5, N auxins = 41, α is set for all concentric circles to be actual circles and η x = η y = 1. The resulting network is depicted in figure 4.11. We depict in figure 4.12 some possible ways of modifying the structure to optimize its performances. It shows even small modifications of the phenotype leads to non-intuitive changes in the genotype that could potentially have a strong effect on heat conduction.

Results and Interest of the method

In the previous paragraph, the method used to optimize the spreader was described. It lies on the combination of two algorithms: the first one, the space colonization algorithm, generate the spreaders while the second one, the genetic algorithm, iteratively finds better solutions to the problem.

Influence of the genetic algorithm parameters

Before discussing the actual results obtained using such algorithm, it needs to be cali- helps us determining whether a certain genetic diversity is kept over generations or if there is a risk of premature convergence to a local optimum. For each population size, the GDA is run three times for statistical purposes and we represent in figure 4.13 the mean value of the fitness among the three runs for each generation and population size. We notice that both maximum and average fitness continuously increase with each generation, which validate the overall genetic algorithm process although increase rate could probably be optimized. Then, comparing the values of maximum thickness, we notice that over time, larger populations lead to higher fitness hence better individuals. This is only related to a better exploration: with more individuals, more alleles can be tested, more mutations can occur hence enabling the algorithm to find a better individual. Quantitatively comparing the evolution of fitness, one can notice a threshold: while all GDA with a population with less than 30 individuals give comparable results, GDA with a population of 50 individuals gives a 10% higher fitness after 50 generations. While the value linked to this threshold is problem dependent and can hardly be compared with data from the literature, it seems that, given the problem complexity, a population of at least 50 individuals is necessary to obtain an acceptable convergence and convergence rate. Comparing now, the average to the maximum fitness it seems that over time, for all population with less than 30 individuals, both are comparable, which may imply that genetic diversity has been lost.

The threshold is even more visible in figure 4.14 in which is represented the evolution of the fitness for different population sizes and generations. A large enough population is necessary to maintain high genetic diversity which enables an acceptable exploration of the optimization space leading to higher fitness. However it should be noted that this is time consuming: for a five individuals population, 2 h were necessary for the GDA to reach fifty generations while it took around 12 h with a population of fifty individuals.

Although it could be seen as a drawback, it seems necessary, here, to use a fifty individuals population.

Fusion front tracking and influence on the heat source temperature

Using the GDA with the parameters described in the previous paragraph, six different structures (one for each population size) are generated after 50 generations and depicted in figure 4.15. Resulting structures are relatively diverse and their respective parameters are reported in Appendix A.3. While in the previous paragraph we compared their fitness, an image of the heat source temperature at the end of the simulation, we want to see here the evolution of the temperature and the influence of both heat spreader structure and fusion front evolution on this temperature. To do so, we represent in figure 4.16 both the die temperature evolution over time and the evolution of the temperature increase rate (temperature time derivative).

In addition to the hybrid structures, we add data obtained with an homogeneous filler having a porosity of 0.8 and notice that the temperature of the heat source for such structure rises way faster. Contrary to the hybrid fillers, heat is directly spread in the effective CHAPTER 4. BEYOND HOMOGENEOUS FILLER heat spreader we notice a four stages as reported in other studies available in the literature [20]. The first stage, only lasts for less than a second and corresponds to the so-called pre-melting stage, when heat is spreading but temperature in the structure remains below the melting temperature (here 108 • C). Then a two steps melting is encountered: in the first step (i.e. for t < 0.5t f visible on figure 4.18) PCM is the vicinity of the heat source is melting. These PCM areas are easily melted due to the relatively high diffusivity of the enhanced PCM. Then, around t < 0.5t f , the fusion front reaches the sides of the optimization domains. At this moment, the fusion front is constrained and PCM further away from the die has to be reached for heat to be stored which causes an increasing the heating rate of the heat source. Finally, around t = 0.8t f , the fusion reaches the end of the domain Ω.

At this point most of the PCM is melted which causes a second rise in the temperature increase rate (seen on figure 4.16b). From this point on, no more heat can be stored through latent heat.

This gives us information on how to improve heat spreader generation for the actual use-case: it should ensure that fusion front remains unconstrained for the longest time possible. Although this is not intuitive, it means that more freedom should be given in the spreader generation without increasing the number of parameters in order to limit the size of the optimization space and time to reach an optimum.

Use-case device design based on Generative design algorithm

In the previous paragraph, we developed and tested an algorithm enabling the growth and optimization of heat spreaders aiming at improving the performance of LHTESS. The GDA was proven efficient optimizing a LHTESS on a simplified 2D domain and is thought to be used on the use-case studied in chapter 3. To do so, it is necessary to adapt the growth and optimization parameter to better accommodate to the new optimization domain.

Spreader generation

As it was done in the previous paragraph, we describe the different parameters defining the growth of the spreader through the space colonization algorithm. First, we need to describe the optimization domain shown in figure 4. 19. It is a 2D projection on the x y plan of the actual LHTESS. As it was done in part 4.2, rather than an actual 3D structure, the spreader is rather obtained through a sweep of a 2D representation which justifies the dimension of the optimization domain.

To properly design the spreader, we use a similar method than the one used in part 4.2: a cylinder of radius R 0 is attached under the die at its center and a given number of branches are connected to it to form the spreader. The design of the branches itself follows a similar principle than in part 4.3.2 but the auxins and source nodes generation is performed using different parameters:

• D: Spatial step size between each node of the network. A higher value of D usually gives more simple or coarser final design, while a smaller value results in a more refined one.

• R 0 : The base cylinder radius to which the branches are attached. It provide a way to rapidly and efficiently spread the heat close to die. If the radius is to large, it may needlessly reduce the space available for the filler.

• N sources : Set the number of source node. Unlike in part 4. • r min : radius of the extremity nodes which is the smallest value found in the network.

Here it is set at 250 µm as it is corresponds to the limit of the additive manufacturing technology.

• r max : radius of the source nodes, that corresponds to the maximum value found in the network.

• C Dist : this parameter picks a single distribution among a list of predefined distributions consisting in concentric half circles oriented in different direction (see fig- • R auxins : corresponds to the radius of the larger concentric half-circle.

• N auxins : set, depending on the auxins distribution (set by C Dist ), either the number of auxins by half-circle or the number of auxins of the larger one. In the last case, the number of auxins on the smaller half-circle is set to keep an even spacing between each auxins on a given half-circle.

• N Circle : Number of concentric half-circles.

• η x and η y : allows to modify the half-circles for them to be half-ellipsoids hence modifying the path for the network growth (0 < η x , η y ≤ 2).

To better illustrate the influence of each parameter, and as it was previously done, we display a reference spreader on figure 4.20 which will be modified afterward. For this arbitrary reference spreader, the following parameters were used: R 0 = 2, N sources = 5, r max = 0.4, C Dist is set for the center of the concentric half-circle centers to be the one In the previous simplified study, it was shown that anytime the fusion front was constrained, it lead to a rise of the temperature increase rate, hence spreader might have to be modified to avoid any constraint. In that sense, additional parameters are added to limit the growth of the spreader to places that are both close to the heat source (i.e. the die) and close to a domain limit. For that two spline curves are defined using two different parameters D 1 and D 2 (see figure 4.22). The growth algorithm is then modified to suppress any auxin nodes that are outside those two curves. In addition, to limit the network growth in some directions to avoid fusion front constraint, it also allows for more conductive material to be found in the lattice filler.

Preliminary 2D study

We have previously shown that for the genetic algorithm to sufficiently explore the optimization domain, given its dimensions, a relatively high number of evaluations is necessary. The same is true considering the optimization of the spreader for our use-case. In the previous paragraph (part 4.3) the GDA was applied to the optimization of a spreader in a 2D optimization domain which limits the time of simulation hence of evaluation of the structure fitness. Going from 2D to 3D actually increases (in our case) the time of simulation by a factor five which practically makes the method hard to apply.

In an attempt to limit simulation time while finding a path towards an optimum heat spreader, it was decided to first perform a preliminary study on the system depicted in to higher increase rate. Actually, once all PCM has melted under the die (t = 0.17 × t f ), a first notable bump is observed before PCM further away from the heat source can be melted. A second notable curve change is observed at t = 0.35 × t f when the fusion front hit the symmetry axis (corresponding to the upper side of the system). The temperature increase rate remains roughly constant until all PCM is melted around t = 0.70 × t f (and above) when it increases again until the end of the simulation. Despite the difference in fitness, we notice that all three optimized spreaders behave in a relatively similar manner, and spreader from the third run gives a difference of around 3 • C at the end of the simulation.

Compared to the spreader free structure, which consists in a simple 2D projection of the samples in part 3 on the x y plan with a porosity of 0.8, we notice a gain of 5 • C which once confirms the interest of topology optimization. If we perform a more thorough analysis of the spreader free structure die temperature profile, we notice an early change of slope that corresponds to the melting of the PCM right underneath the die. As no spreader part is found underneath the die in this case, more PCM can be melted in a shorter time.

However, once all PCM has been melted underneath, heat spreads in the enhanced PCM without any help of a spreader which causes a higher temperature increase rate after this first step. Comparing now the evolution of the fusion front, we see that we are able to direct it in some directions of space to avoid some constraints. It is particularly visible between t = 0.35 × t f and t = 0.50 × t f where temperature increase rate is rising up in the spreader free case while it is maintained almost constant for the spreader based devices. Comparing the three optimized spreaders, very similar behaviors are observed. This might be due to the relative lack of topology diversity which combined with the high input power that may hide any difference, at least during phase transition. Comparing them to the non-optimized spreader, a similar trend appears. This could signify two things :

either in these conditions of high input power on a relatively short time, a refined heat spreader topology is of no help and the only important parameter is the amount of highly conductive material (here aluminum) in close direct contact with the die or the GDA could not find a global optimum given the set of parameters used here. As shown in part 4. Material with Penalization method) has already been proven efficient in a variety of cases [6; 10; 22]. If no optimal solution can be found through SIMP, this signifies that the theory of non-optimality of tree-like structures over laminar structure (like the non-optimized spreader in part 4.2) may actually be verified; in this case a more conventional parametric analysis on a restricted number of parameters may be the best solution. All these remain open question to this day.

Conclusion

After finding out, in the second and third chapters, the most suitable simulation assumptions and modeling tools, this fourth chapter aimed at optimizing the use-case device in order to reach the target die temperature of 120 • C while maintaining the same level of mass reduction.

In the first part of the chapter, following the cascaded storage theory, two different heat spreading structures made of aluminum are added to the PCM enclosure of the usecase device. To avoid increasing the total mass of the device and reducing the amount

Conclusion & Perspectives

The use of PCM is thought to be a promising solution for passive cooling of power electronics as long as efficient thermal conductivity enhancement is provided. In this study, enhancement is based on the use of additively manufactured structure which topology can be tuned to better adapt to a variety of cases. The advantage of such method relies in its ability to adapt the filler to PCM ratio locally in order to optimize heat transfer and PCM melting rate to limit as much as possible the dies temperature rise.

After a quick presentation of the PCM global properties and their diversity, chapter 1 was dedicated to the study of state-of-the-art PCM modeling and enhancement which drew the baseline of the entire study. It introduced the enthalpy method that allows to model fusion without precise description of the solid/liquid interface, reducing simulation computational cost. Finally, enhancement techniques based on the use of porous conductive materials as well as way to simplify the description of the resulting composite material were presented.

Based on the study of the current thermal conductivity enhancers (i.e. metal foam) we proposed, in chapter 2, a model of effective thermal conductivity for lattice structures.

Actually, studies on foams showed that more homogeneous and regular structures could increase the thermal conductivity of the filler due to the elimination of a variety of thermal resistances. The development of the model revealed notable progress as thermal conductivity was increased by up-to 75% in some directions of space compared to foam. While this first part of the study was performed on an isolated cell, the second part concentrated on the integration of such dissipating structures in actual devices by considering side effects. Actually, those are additively manufactured which comes with limitation regarding the minimal dimensions of the different parts. To maintain a sufficiently high level of porosity, a reduced number of stacked cells might arise which compromises homogenization. It was shown that these are due to constriction phenomena, which arise when heat flows between two solids of different section area. In this situation, flow is constrained (or widen) which induces additional thermal resistances, reducing the actual Effective Thermal Conductivity (ETC). Constriction was implemented into the above mentioned model and allowed a more refined ETC estimation. Using this improved model, additional pa-

CONCLUSION & PERSPECTIVES

rameters of the porous network are accounted for, like for example pore size, number of stacked cells or heat spreader dimensions. Both simplified and constriction models were finally validated experimentally, using flash laser method which showed good agreement.

Using the models and structures previously developed, chapter 3 concentrated on the optimization of lattice structure infused with PCM based device for an actual use-case in terms of dimensions, aluminum to PCM volume ratio and PCM choice. The first step in the optimization process consisted in characterizing a series of PCM in terms of melting temperature, latent heat of fusion and specific heat profile. This characteristic were then used in the enthalpy model to accurately simulate fusion. Employing the different models, a parametric study was run to determine the optimum material distribution in the use-case, which showed that maximizing the quantity of PCM in the device could be detrimental to its performance. Actually, it was demonstrated that an optimum PCM quantity existed that allowed efficient heat spreading close to the heat source, heat that could be then as efficiently stored in the PCM further away. Such optimization helped develop devices that could potentially decrease end-of-mission temperature as well as device mass by up to 20% compared to existing solutions. In the last part of the chapter, a series of such devices were manufactured to confirm simulation assumptions. Although the experiments confirmed the trend already depicted by the simulations, they showed that the assumptions were over-simplified and that refined ones (i.e. LTNE) had to be considered.

This assumption led to higher temperature than expected, demonstrating the need for more advanced solutions.

Given the results of the previous chapter, chapter 4 proposed a method that helped further improve heat spreading and storing at the same time : Cascaded storage. This method relies on a smart way of distributed both conductive (i.e. aluminum) and storing (i.e. PCM) in the device in order to maximize heat spreading and homogenize the temperature. Here, this was performed by introducing a so-called heat spreader inside the improved PCM enclosure to better spread heat and avoid hot points close the heat source.

To maximize the benefit of such technique the shape of the spreader was optimized by the combination of two algorithms: space colonization and genetic. The former, given a set of parameters, generates bio-inspired structure naturally adapted to heat conduction while the latter optimizes the parameter set in order to minimize the heat source temperature. Although promising results were obtained in 2D, the algorithm failed to improve basic solutions in 3D which might be due to limited set of parameters used to generate the structure. Even though more parameters could be considered, it could be detrimental in terms of calculation time and different optimization methods need to be assessed.

This work has confirmed the potential of PCM based devices as efficient solution for the thermal management of power electronics. Using advanced filler, additively manufac-

CONCLUSION & PERSPECTIVES

tured, in combination with Cascaded storage principles, a temperature gain of 14% and a mass gain of 20% were observed both theoretically and experimentally. Although promising, some points need to be further studied to improve both performances of the device and understanding of the physical processes involved in the spreading and storing. Those points are listed below:

• Values of integral heat transfer coefficient in chapters 3 and 4 had to be estimated by simulation to experiments correlation which may lead to error as they could include another physical phenomenon interpreted as thermal lagging. One way to overcome this issue would be to extend the closure problem described by Quintard et al. to a variety of topologies to improve estimation;

• As mentioned previously, generative design algorithm may not be adapted to the optimization of such structure. Actually, the number of parameters describing it might be either to low to fully optimize its design or to high for the genetic algorithm combined with an FEA model to be efficient enough. Different methods of generation could be used in order to improve the optimization process: like Solid-Isotropic Structure with Penalization (SIMP) which relies on the definition of a vector defining the given ratio of conductive material or PCM for individual voxels (or pixels in 2D) on a given domain.

• Progress in additive manufacturing can help improving device performances in two ways. The first one is related to the resolution: scaling down the minimum strut thickness could help reduce the pore size for a given porosity which increases the specific area, lessening thermal lagging which has been shown to drastically impact device performance and to a lesser extent thermal constriction as mentioned in chapter 2. Progress in terms of printed material could also help improving the device : for now, only some aluminum alloys can be used in additive manufacturing (mainly AS7) which thermal properties are usually less than those of Al6061, for example.

A.1.4 PCM Characterization

4. 1 A. 1

 11 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Preliminary study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 2D Generative design algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Use-case device design based on Generative design algorithm . . . . . . . . 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Supplementary figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I 1 Phase Change Material (PCM) melting in an idealized metal foam filler. . . . 1.1 PCM Classification [22]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Melting temperatures and Latent Heat of Fusion of existing PCM [10]. . . . . 1.3 PCM main applications: (a) Thermal Damping, (b) Limiting Thermal Excursion, (c) Maintaining a constant temperature. . . . . . . . . . . . . . . . . . . 1.4 Photography of the Dover House designed by Eleanor Raymond [29]. . . . . 1.5 Heating system of the Dover House designed by Maria Telkes [29]. . . . . . . 1.6 Evolution of the temperature distribution at five different instants t 1 < ... < t 5 using the one-dimensional definition of the Stefan Problem with Dirichlet condition in x = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Evolution of the phases temperature difference as a function of : (a) . . . . . 1.8 Computational domain and representative mesh for Direct Numerical Simulation (DNS) [73]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Evolution of the melted fraction as a function of F o for the four different cases simulated by Yang et al. [73]. . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Comparison between solid-liquid interface evolution obtained separately with volume average method and DNS: (a) Considering fusion-induced convection, (b) Not considering fusion-induced convection [73]. . . . . . . . . . 1.11 Comparison between temperature distribution evolution obtained separately with volume average method and DNS: (a) Considering fusion-induced convection, (b) Not considering fusion-induced convection [73]. . . . . . . . . . 1.12 PCM Melting front comparison at four selected times under inclination angle of 0°, 30°, 60°, 90°[75]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 PCM Melting fraction comparison at four selected times under inclination angle of 0°, 30°, 60°, 90°[75]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.14 Comparison of specific Effective Thermal Conductivity (ETC) for different types of filler [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.15 Examples of Architected Enhancers: (a) Honeycomb, (b), Truss-lattice structure, (c), TPMS [89] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii LIST OF FIGURES 1.16 Foam-like and TPMS structures studied by Ahmed et al. [86]: (a) Tetrakaidecahedron, (b) Primitive, (c) IWP and (d) Gyroid. . . . . . . . . . . . . . . . . . . 1.17 Parallel (on the left) and Series (on the right) models comparison. . . . . . . 1.18 Scanning Electronic Microscope image of an aluminum foam displaying lumping at the node [7]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.19 Calmidi and Mahajan [95] foam representative medium: (a) Hexagonal Network and (b) Unit Cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.20 Bhattacharya representative medium unit cell [43]. . . . . . . . . . . . . . . . 1.21 Foam three-dimensional representative network unit-cell according to: (a)

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.22 Boomsma and Poulikakos unit cell : (a) Tetrakaidecahedron with cylindrical ligaments and cubic nodes, (b) Geometrical breakdown of the unit cell of the Tetrakaidecahedron. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.23 Yang et al. Foam Unit-cell [7]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.24 Foam ETC analytical models comparison. . . . . . . . . . . . . . . . . . . . . . 1.25 Strut profile and cross section variation of (a) a Urethane foam and (b) an Aluminum foam obtained for X-ray tomography [100]. . . . . . . . . . . . . . 1.26 (a) Voronoi tessellation generated from seed positions produced by an arrangement of random seed points (b) Laguerre Voronoi Tessellation generated using the same arrangement of random seed points [101]. . . . . . . . . 1.27 Normalized effective thermal conductivity of random monomodal cells (black triangles with dash curve) and the associated porosity of the Control Volume (CV) (open circles with dot curve) versus ratio of strut cross-sectional areas [104]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.28 Influence of the pore size distribution on the network . . . . . . . . . . . . . 1.29 Zheng et al. optimum porosity gradient arrangement. From left to right: Uniform porosity, Horizontal gradient, Vertical gradient, Two dimensional gradient [107]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.30 Pin-foam enhanced PCM system [12]. . . . . . . . . . . . . . . . . . . . . . . . 2.1 Lattice cells and dimensions : (a) BCC, (b) FCC, (c) BCC z , (d) FCC z . . . . . . . . 2.2 Lattice Unit Cell Slicing (BCC). . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Lattice Unit Cell Slicing (FCC): (a) out-of-plane direction, (b) in-plane direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Lattice Unit Cell Slicing (BCC z ): a out-of-plane direction, b in-plane direction. 2.5 Lattice Unit Cell Slicing (FCC z ): a out-of-plane direction, b in-plane direction. 2.6 Identified node types [10]: a BCC node, b BCC z node, c FCC node. . . . . . . . LIST OF FIGURES 3.2 Influence of the sample mass on the resulting heat flow measurement for two PureTemp 108 samples of respective mass m 1 = 3.83 mg and m 2 = 8.52 mg.127 3.3 Resulting specific heat capacity profiles using a Differential Scanning Calorimetry (DSC) dynamic method on a PureTemp 108. . . . . . . . . . . . . . . . . . 3.4 DSC Isothermal step method realized on a PureTemp 108 sample : (a) Temperature program and resulting heat flow, (b) Stored energy at the different temperature steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 PCM specific heat profile derived from DSC experiment: (a) Model fit and comparison with the experimental data for the PureTemp 108, (b) : Comparison of the different PCM resulting profiles. . . . . . . . . . . . . . . . . . 3.6 Use-case PCM based system definition: (a) Isometric view, (b) Section view. 3.7 Geometry parametric analysis results : (a) End-of-mission die temperature and (b)End-of-mission PCM Melting ratio. . . . . . . . . . . . . . . . . . . . . . 3.8 End-of-mission die temperature for varying side walls and top cover thickness.137 3.9 Influence of PCM choice on the dies temperature profile. . . . . . . . . . . . . 3.10 Geometry parametric analysis results using Rubitherm RT 100HC instead of Rubitherm RT 90HC : (a) End-of-mission die temperature and (b)End-ofmission PCM Melting ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 View of sample #3 using tomography images: (a) View of the top of the sample, (b) View underneath the top cover (sealing plug baseplates are seen on the sample edges) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Experimental setup epoxy board designs: (a) Bottom layer, (b) Middle layer and (c) Top layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Model validation experimental setup. . . . . . . . . . . . . . . . . . . . . . . . 3.14 Bare aluminum sample model used in simulations. . . . . . . . . . . . . . . . 3.15 Experimental setup Calibration : Comparison of the predicted and measured die temperature profiles for bare aluminum samples of different height (i.e. 9, 10 and 11 mm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.16 Simulation to experiment correlations for varying pore size and local thermal equilibrium assumption: (a) Temperature profiles and (b) Relative error compared to experimental data. . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 17

 17 Simulation to experiment correlations for varying porosity and local thermal equilibrium assumption: (a) Temperature profiles and (b) Relative error compared to experimental data. . . . . . . . . . . . . . . . . . . . . . . . . . .3.18 Simulation to experiment correlations for varying overall package thickness and local thermal equilibrium assumption: (a) Temperature profiles and (b) Relative error compared to experimental data. . . . . . . . . . . . . . . . . . . xvi LIST OF FIGURES 3.19 Simulation to experiment correlations for varying top cover thickness and local thermal equilibrium assumption: (a) Temperature profiles and (b) Relative error compared to experimental data. . . . . . . . . . . . . . . . . . . . . 3.20 Assessment of the Local Thermal Non-Equilibrium (LTNE) assumption on the die temperature profile considering use-case boundary conditions. . . . 4.1 Feng et al. tested fillers [1] : (a) Metallic foam, (b) Array of fins and (c) Finned foam. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Hybrid Lattice Filler devices top view (Enhanced PCM in blue, Aluminum in grey) : (a) Fin and bare PCM and (b) Fin and adapted lattice filler. . . . . . . . 4.3 Hybrid filler comparison with existing solutions using: (a) Local Thermal Equilibrium (LTE) assumption (average and maximum die temperature) or (b) LTNE assumption (average and maximum die temperature). . . . . . . . . 4.4 Preliminary study optimization domain. . . . . . . . . . . . . . . . . . . . . . 4.5 Preliminary study domain filled with a heat spreader and enhanced PCM. . . 4.6 Space colonization algorithm principle. . . . . . . . . . . . . . . . . . . . . . . 4.7 Thickening algorithm principle. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Alpha shape creation process [16]. . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Genetic algorithm flow chart. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Comparison of heuristic and genetic algorithms evolution. . . . . . . . . . . 4.11 Reference parameter resulting (a) Network and (b) Heat spreader. . . . . . . 4.12 Evolution of the resulting network from the reference one when: (a) D = L X /240, (b) N source = 13, (c) N auxins = 21, (d) r max = L x /60, (e) Half-circles undergo a 180°rotation, (f) N circle = 8, (g) α is modified to extend the colonization space, (h) η x = 0.6 and (i) η y = 0.6. . . . . . . . . . . . . . . . . . . . . 4.13 Resulting (a) maximum and (b) average fitness for different population sizes over 50 generations running the Generative Design Algorithm (GDA). . . . . 4.14 Resulting maximum fitness for different population sizes. . . . . . . . . . . . 4.15 Heat spreaders generated by GDA after 50 generations using population sizes of: (a) 5 individuals, (b) 10 individuals, (c) 15 individuals, (d) 20 individuals, (e) 30 individuals and (f) 50 individuals. . . . . . . . . . . . . . . . . . . . . . . 4.16 Die temperature profile (a) and Heating rate (b) for the six GDA optimized spreaders. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 GDA resulting structure compared temperature distribution and average temperature (T avg ) evolution over time. . . . . . . . . . . . . . . . . . . . . . . . . 4.18 GDA resulting structure compared melting ratio evolution over time. . . . . . 4.19 Use-case optimization space for GDA. . . . . . . . . . . . . . . . . . . . . . . . 4.20 Reference heat spreader generated in the use-case optimization domain. . . xvii LIST OF FIGURES 4.21 Evolution of the resulting network from the reference one when: (a) R 0 = 5, (b) N source = 8, (c) Half-circles undergo a 180°rotation, (d) R Source = 2R 0 , (e) N auxins = 20 , (f) N circle = 10, (g) η x = 1.75 and (h) η y = 1.75 . . . . . . . . . . . 182 4.22 Spreader growth adaptation to avoid fusion front constraint. . . . . . . . . . 183 4.23 GDA (a) maximum and (b) average resulting finesses for the use-case 2D preliminary study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 4.24 (a) Temperature and (b) Temperature increase rate profiles for the three optimized spreader and spreader free structure. . . . . . . . . . . . . . . . . . . 185 4.25 Optimized spreaders 2D fusion front tracking over time. . . . . . . . . . . . . 186 4.26 Compared top cover influence on the end-of-mission temperature for spreader free or optimized spreader based devices. . . . . . . . . . . . . . . . . . . . . . 187 4.27 Compared average die temperature profiles of devices using hybrid (optimized and non-optimized) and lattice fillers. . . . . . . . . . . . . . . . . . . . 188 A.1 Model of ETC comparison and theoretical validation with : (a) a BCC cell, (b) a BCC z cell, and (c) a FCC z cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . I A.2 Calorimeter calibration using a POCO sample (courtesy of Influtherm). . . . VI A.3 Influence of the sample mass on the resulting heat flow measurement for two Rubitherm RT80HC samples of respective mass m 1 = 4.02 mg and m 2 = 9.18 mg. VII A.4 Influence of the sample mass on the resulting heat flow measurement for two Rubitherm RT90HC samples of respective mass m 1 = 3.83 mg and m 2 = 8.52 mg. VII A.5 Influence of the sample mass on the resulting heat flow measurement for two Rubitherm RT100HC samples of respective mass m 1 = 3.83 mg and m 2 = 8.52 mg.VIII A.6 DSC Isothermal step method realized on a Rubitherm RT80HC sample: (a) Temperature program and resulting heat flow, (b) Stored energy at the different temperature steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VIII A.7 DSC Isothermal step method realized on a Rubitherm RT90HC sample: (a) Temperature program and resulting heat flow, (b) Stored energy at the different temperature steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IX A.8 DSC Isothermal step method realized on a Rubitherm RT100HC sample: (a) Temperature program and resulting heat flow, (b) Stored energy at the different temperature steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IX xviii List of Tables

CHAPTER 1 .

 1 Figure 1.1: PCM Classification [22].

Figure 1 . 2 :

 12 Figure 1.2: Melting temperatures and Latent Heat of Fusion of existing PCM [10].

Figure 1 . 3 :

 13 Figure 1.3: PCM main applications: (a) Thermal Damping, (b) Limiting Thermal Excursion, (c) Maintaining a constant temperature.

CHAPTER 1 .Figure 1 . 5 :

 115 Figure 1.5: Heating system of the Dover House designed by Maria Telkes [29].

Figure 1 . 7 :

 17 Figure 1.7: Evolution of the phases temperature difference as a function of : (a) their respective thermal conductivities ratio, (b) the fluid volume ratio.

Figure 1 . 9 :

 19 Figure 1.9: Evolution of the melted fraction as a function of F o for the four different cases simulated by Yang et al. [73].

Figure 1 . 10 :

 110 Figure 1.10: Comparison between solid-liquid interface evolution obtained separately with volume average method and DNS: (a) Considering fusion-induced convection, (b) Not considering fusion-induced convection [73].

Figure 1 . 11 :Figure 1 . 12 :

 111112 Figure 1.11: Comparison between temperature distribution evolution obtained separately with volume average method and DNS: (a) Considering fusion-induced convection, (b) Not considering fusion-induced convection [73].

igure 1 . 13 :

 113 PCM Melting fraction comparison at four selected times under inclination angle of 0°, 30°, 60°, 90°[START_REF] Yang | Effect of inclination on the thermal response of composite phase change materials for thermal energy storage[END_REF].

Figure 1 . 14 :

 114 Figure 1.14: Comparison of specific ETC for different types of filler [5].

Figure 1 . 15 :

 115 Figure 1.15: Examples of Architected Enhancers: (a) Honeycomb, (b), Truss-lattice structure, (c),

  et al. [86] compared the efficiency of foam-like enhancer in the form of a combination of tetrakaidecahedron (figure 1.16a) to different types of architected materials based on Triply Periodic Minimal Surface (TPMS). TPMS are non self-intersecting surface in R 3 , which has a crystalline structure repeated in three independent directions and has zero mean curvature at each point [91]. To assess the performances of those structures Ahmed et al. first measured their ETC and the one of the foam-like structure for porosity of 90%. It was determined that the Primitive (figure 1.16b) structure had the highest ETC with 11.70 W•K -1 •m -1 then came the (figure 1.16c) with 11.16 W•K -1 •m -1 , the Gyroid (figure 1.16d) with with 10.46 W•K -1 •m -1 and finally the foam-like structure with 7.09 W•K -1 •m -1 which shows the advantage of such structure. To further deepen the understanding of such structure behavior, Ahmed et al. ran dynamic simulation in which the melting of the PCM is monitored. Results showed that foam-like structure was still

Figure 1 .Figure 1 . 17 :

 1117 Figure 1.16: Foam-like and TPMS structures studied by Ahmed et al. [86]: (a) Tetrakaidecahedron, (b) Primitive, (c) IWP and (d) Gyroid.

Figure 1 . 18 :Figure 1 . 19 :

 118119 Figure 1.18: Scanning Electronic Microscope image of an aluminum foam displaying lumping at the node [7].

Figure 1 . 21 :Figure 1 . 22 :

 121122 Figure 1.21: Foam three-dimensional representative network unit-cell according to: (a) LordKelvin (Tetrakaidecahedron)[START_REF] Thomson | On the Division of Space with Minimum Partitional Area[END_REF], (b) Weaire and Phelan[START_REF] Weaire | A counter-example to kelvin's conjecture on minimal surfaces[END_REF] and (c) Gabbrielli[START_REF] Gabbrielli | A new counter-example to Kelvin's conjecture on minimal surfaces[END_REF].

  first used the Tetrakaidecahedron as a unit cell (see figure 1.22a) for the determination of the foam ETC. For simplicity purpose, it is represented with cylindrical ligament and cubic nodes. Boomsma and Poulikakos used a similar method than Calmidi and Mahajan with a 3D unit-cell depicted on figures 1.22b that is sliced into four different layers on which the Parallel model is applied. Boomsma and Poulikakos first expressed the solid volume in each layer as a function of both d = a/L and e = r /L, leading to:

CHAPTER 1 .Figure 1 . 23 :

 1123 Figure 1.23: Yang et al. Foam Unit-cell [7].

2 . ( 1 . 90 )

 2190 To calculate the integral, Yang et al. sliced the cell into three parts: one from the bottom of the cell to the top of the nodes at the bottom (z ∈ [0, t /2]), the second from the top of the node at the bottom to the bottom of the node at the top of the cell (z ∈ [t /2, 2L/2t /2]) and finally from the bottom of the top nodes to the top of the cell(z ∈ [ 2L/2t , 2L/2]

. 93 )

 93 Substituting equation 1.93 into equation 1.88, it gives:

  with a = 327.25811, b = -1075.55645, c = 1182.83207 and d = -434.55535. Yang et al. latter compared the value of ETC obtained using Dai et al. improved model and their own model to experimental data and showed that while Dai et al. model exhibits a 13.2% RMS deviation, their own model only exhibits a 5.3% RMS deviation. On figure 1.24, the previous models are compared to experimental values to assess their respective precision. It is observed that Parallel and Series models are respectively an overestimate and an underestimate of the foam ETC. Then, comparing the other model in terms of RMS deviation, the different model can be ranked as follows:

Figure 1 . 25 :

 125 Figure 1.25: Strut profile and cross section variation of (a) a Urethane foam and (b) an Aluminum foam obtained for X-ray tomography [100].

  [START_REF] Jang | On the microstructure of open-cell foams and its effect on elastic properties[END_REF]. Jang et al. studied foams made out of different materials (mainly polyester urethane and Duocel aluminum) with different porosity and average cell size using X-ray tomography to establish statistically the cell size, ligament length distribution, material distribution along the ligament, geometry of the node and cell anisotropy. From this study Jang et al. determined that both cell size and ligament length inside the network follow a log-normal distribution with a relatively tight peak for both urethane and aluminum foams. Then, Jang et al. studied the material distribution along the foam strut and determined that during the foam process, material tend to concentrate at the node to minimize their specific area forming hyperbolic strut rather than cylindrical (see figure1.25). The cross section shape is different depending on the material and foaming process, as a consequence it is observed that the urethane foam has a three-cusps hypocycloid cross section of Plateau borders (figure 1.25a) while the aluminum foam's struts have a rounder profile (figure 1.25b).

Figure 1 .

 1 Figure 1.26: (a) Voronoi tessellation generated from seed positions produced by an arrangement of random seed points (b) Laguerre Voronoi Tessellation generated using the same arrangement of random seed points [101].

Figure 1 .Figure 1 . 27 :

 1127 Figure 1.26 shows how a 2D domain can be divided using such method: in figure 1.26a, a basic Voronoï tessellation is used leading to a relatively homogeneity in the region size distribution while on figure 1.26b a weighted version of this algorithm is used showing a more poly-dispersed division. Using this technique, a number of studies [45; 104; 105] were realized to determined the influence of the different parameters defined by Jang et al.. The first parameter to be studied was the porosity that varied in the studied papers from 0.75 to 0.96. It was shown that similarly to what was predicted by analytical model, the ETC does not vary linearlywith the porosity and still fit experimental data proving once again the efficiency of such analytical models. The second parameter to be studied was the strut profile and the concentration of matter at the nodes. It was shown that as the node to strut diameter ratio increased, the global ETC of the foam was reduced which can be explained by the additional thermal resistance generated in the middle of the strut where the cross sectional area is the smallest. This is visible on figure 1.27, on which is represented the ETC ratio between hyperbolic strut foam with varying curvature and cylindrical strut foam depend-

Figure 1 . 28 :

 128 Figure 1.28: Influence of the pore size distribution on the network ETC [45].

Figure 1 . 29 :

 129 Figure 1.29: Zheng et al. optimum porosity gradient arrangement. From left to right: Uniform porosity, Horizontal gradient, Vertical gradient, Two dimensional gradient [107].

  [12] studied a system in which a foam of various porosity, material and porosity gradient (in the direction of the fins) is combined with an array of fins oriented vertically. The system depicted in figure1.30 is subjected to a constant temperature at its bottom, natural convection on its top and is isolated on the sides. Initially, the PCM is already melted and Yang et al. observed its solidification. Unlike in Zheng et al. study, it should be noted that the problem is conduction dominated. Yang et al. observed similar results than Zheng et al., as the introduction of the pin array into the metal foam PCM increased the solidification rate by 30.5% and up to 34.7% for pin-foam hybrid structure with a graded porosity (with the lowest porosity at the bottom of the structure).

  Figure 2.1: Lattice cells and dimensions : (a) BCC, (b) FCC, (c) BCC z , (d) FCC z .

CHAPTER 2 .Figure 2 . 2 :

 222 Figure 2.2: Lattice Unit Cell Slicing (BCC).

. 21 )Figure 2 . 3 :

 2123 Figure 2.3: Lattice Unit Cell Slicing (FCC): (a) out-of-plane direction, (b) in-plane direction.

. 28 )Figure 2 . 4 :

 2824 Figure 2.4: Lattice Unit Cell Slicing (BCC z ): a out-of-plane direction, b in-plane direction.

Figure 2 . 5 :

 25 Figure 2.5: Lattice Unit Cell Slicing (FCC z ): a out-of-plane direction, b in-plane direction.

  FCC z Cell: Porosity and Topology Parameters Calculation

Figure 2 . 6 :

 26 Figure 2.6: Identified node types [10]: a BCC node, b BCC z node, c FCC node.

. 60 )

 60 Comparing Vaissier et al. work to the the structures defined in figure 2.1 we identify three types of nodes based on their so-called valency (number of strut linked to them) CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING and strut agencement around them:

Figure 2 . 8 :

 28 Figure 2.8: Normalized ETC validation: (a) k z , (b) k x y .

CHAPTER 2 .Figure 2 . 9 :

 229 Figure 2.9: Strut radius influence on porosity.

Figure 2 . 10 :

 210 Figure 2.10: Topology parameters comparison:(a) G s,z , (b) G s,x y .

Yang's model

  is a best case scenario. As a matter of fact, Yang et al. considered a perfect tetrakaidecahedron cell to develop their model for simplicity purpose though it has be shown by Baillis et al.[3] that foam random features reduce a real foam ETC. This reinforces lattice structures advantages over foam as PCM enhancers.The choice of one topology over the others may depend on several parameters : geometry of the system in relation with position of the heat source, manufacturing issues, thermal contact resistance between the composite and other part of the thermal management device (see part 2.3.1). So, even if G s,z and G s,x y are of crucial importance regarding the choice of cell for a given application, a last parameter needs to be accounted for. Actually, CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

Figure 2 . 11 :

 211 Figure 2.11: Additional constriction resistance at the perfect interface between two cylinders of very different cross sectional area.

Figure 2 . 12 :

 212 Figure 2.12: Observation of thermal constriction at the interface between the enhanced PCM and the heat spreaders.

Figure 2 . 14 :

 214 Figure 2.14: Gladwell's micro-contact modeling: (left) Isometric view, (right) Top view.

  i ) i ∈N : Constants allowing a variety of Neumann or Dirichlet boundary conditions, -Γ, Λ: Functions depending on the boundary conditions, -C n , D n : Series depending on the boundary conditions, -J 0 : Bessel function of the first kind and first order, -I 0 : Modified Bessel function of the first kind and first order, -ξ, λ n : Integration variables (m -1 ).

CHAPTER 2 .Figure 2 . 15 :

 2215 Figure 2.15: Dimensionless constriction resistance in relation with parameters β and η.

Figure 2 . 18 :

 218 Figure 2.18: Model of ETC comparison and theoretical validation with a FCC cell.

Figure 2 .

 2 Figure 2.8 already depicted the influence of both cell topology and porosity on ETC.

CHAPTER 2 .

 2 PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

Figure 2 . 19 :

 219 Figure 2.19: Influence of the combined effect of strut radius and heat spreader thickness.

CHAPTER 2 .Figure 2 .

 22 Figure 2.20: Sample characteristic dimensions.

CHAPTER 2 .Figure 2 . 21 :

 2221 Figure 2.21: Main SLM parameters [52].

ure 2 .Figure 2 . 22 :

 2222 Figure 2.22: Additively manufactured samples defects: (a) Roughness observed on a type 3 sample, (b) Holes observed on the heat spreader of a type 7 sample.

CHAPTER 2 .Figure 2 . 24 :

 2224 Figure 2.24: Base Material Coefficient of Thermal Expansion Measurement.

CHAPTER 2 . 1 )- 5 Figure 2 . 25 :

 215225 Figure 2.25: Base Material Relative Elongation and Coefficient of Thermal Expansion as functions of the Temperature.

Figure 2 . 26 :

 226 Figure 2.26: Flash Laser Method Machine: (a) Global View, (b) Sample Chamber View.

CHAPTER 2 .Figure 2 . 27 :

 2227 Figure 2.27: Different types of thermocouples scheme: (a) Hot-Junction Thermocouple, (b) Hotjunction-less Thermocouple.

CHAPTER 2 .Figure 2 . 28 :

 2228 Figure 2.28: Experimental and modeled thermogram obtained for the base material sample (courtesy of Influtherm).

Figure 2 . 29 :

 229 Figure 2.29: Model comparison for flash laser method thermogram fit: (a) Sample 3 fit using finite difference model, (b) Sample 3 fit using Parker model,(c) Sample 5 fit using finite difference model, (d) Sample 5 fit using Parker model.

Figure 2 . 30 :

 230 Figure 2.30: Sample 3 five layers model.

  .114) with x the position defined in figure 2.30, t the time and h r a convective exchange coefficient with the outside environment that is at temperature T ext . On top of the sample (x = 0), we define the following boundary condition (B.C.) accounting for the flash applied on the top surface:

1 3 . 4

 134 Differential Scanning Calorimetry Dynamic Method . . . . . . . . . 125 3.2.2 Differential Scanning Calorimetry Isothermal Step Method . . . . . 128 3.2.3 DSC Methods Combination and Effective Specific Heat Capacity Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 3.3 Use-case Parametric Study . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 3.3.1 Use-case Definition and Physical Model . . . . . . . . . . . . . . . . 132 3.3.2 Geometry Parametric Analysis . . . . . . . . . . . . . . . . . . . . . . 134 3.3.3 PCM Choice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 3.3.4 Comparison with different types of filler . . . . . . . . . . . . . . . . 140 3.3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 Experimental Validation of the Model Implementation . . . . . . . . . . 141 3.4.1 Samples description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 3.4.2 Experimental Setup Description . . . . . . . . . . . . . . . . . . . . . 144 3.4.3 Experimental Setup Calibration . . . . . . . . . . . . . . . . . . . . . 146 3.4.4 Experimental Results and Discussion . . . . . . . . . . . . . . . . . . 148 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 3.6 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

- 1 Figure 3 . 1 :

 131 Figure 3.1: DSC Programs for different heating ramps (i.e. 1 K•min -1 , 5 K•min -1 and 10 K•min -1 ), realized on a PureTemp 108 sample of mass 3.83 mg and resulting Heat Flow measurement.

Figure 3 . 2 :

 32 Figure 3.2: Influence of the sample mass on the resulting heat flow measurement for two PureTemp 108 samples of respective mass m 1 = 3.83 mg and m 2 = 8.52 mg.

CHAPTER 3 .- 1 .K - 1 ) 4 H- 1 .K - 1 ) 4 Figure 3 . 3 :

 311411433 Figure 3.3: Resulting specific heat capacity profiles using a DSC dynamic method on a PureTemp 108.

Figure 3 . 4 :

 34 Figure 3.4: DSC Isothermal step method realized on a PureTemp 108 sample : (a) Temperature program and resulting heat flow, (b) Stored energy at the different temperature steps.

- 1 .K - 1 ) 4 HrFigure 3 . 5 :

 11435 Figure 3.5: PCM specific heat profile derived from DSC experiment: (a) Model fit and comparison with the experimental data for the PureTemp 108, (b) : Comparison of the different PCM resulting profiles.

CHAPTER 3 .Figure 3 . 6 :

 336 Figure 3.6: Use-case PCM based system definition: (a) Isometric view, (b) Section view.

Figure 3 . 7 :

 37 Figure 3.7: Geometry parametric analysis results : (a) End-of-mission die temperature and (b)Endof-mission PCM Melting ratio.

Figure 3 . 8 :

 38 Figure 3.8: End-of-mission die temperature for varying side walls and top cover thickness.

Figure 3 . 9 :

 39 Figure 3.9: Influence of PCM choice on the dies temperature profile.

Figure 3 . 10 :

 310 Figure 3.10: Geometry parametric analysis results using Rubitherm RT 100HC instead of Rubitherm RT 90HC : (a) End-of-mission die temperature and (b)End-of-mission PCM Melting ratio.

CHAPTER 3 .Figure 3 . 11 :

 3311 Figure 3.11: View of sample #3 using tomography images: (a) View of the top of the sample, (b) View underneath the top cover (sealing plug baseplates are seen on the sample edges) .

CHAPTER 3 .Figure 3 . 12 :

 3312 Figure 3.12: Experimental setup epoxy board designs: (a) Bottom layer, (b) Middle layer and (c) Top layer.

Figure 3 . 13 :

 313 Figure 3.13: Model validation experimental setup.

Figure 3 . 14 :

 314 Figure 3.14: Bare aluminum sample model used in simulations.

HFigure 3 . 15 :

 315 Figure 3.15: Experimental setup Calibration : Comparison of the predicted and measured die temperature profiles for bare aluminum samples of different height (i.e. 9, 10 and 11 mm).

HHFigure 3 . 16 :

 316 Figure 3.16: Simulation to experiment correlations for varying pore size and local thermal equilibrium assumption: (a) Temperature profiles and (b) Relative error compared to experimental data.

HHFigure 3 . 17 :

 317 Figure 3.17: Simulation to experiment correlations for varying porosity and local thermal equilibrium assumption: (a) Temperature profiles and (b) Relative error compared to experimental data.

HHFigure 3 . 18 :

 318 Figure 3.18: Simulation to experiment correlations for varying overall package thickness and local thermal equilibrium assumption: (a) Temperature profiles and (b) Relative error compared to experimental data.

HHFigure 3 . 19 :

 319 Figure 3.19: Simulation to experiment correlations for varying top cover thickness and local thermal equilibrium assumption: (a) Temperature profiles and (b) Relative error compared to experimental data.

Figure 3 . 20 :

 320 Figure 3.20: Assessment of the LTNE assumption on the die temperature profile considering usecase boundary conditions.

  their work [1], Feng et al. discuss the possible use of different fillers for the enhancement of PCM intended to be used as energy storage solution. Three fillers were used in this study: a foam, an array of fins and an hybrid filler made of a combination of both foam and array of fins (see figure 4.1). Feng et al. tested the different filler on a similar application and found out that the hybrid filler outperformed the other ones. To compare the different filler, Feng et al. calculates the heat transfer coefficient h defined as follow:

Figure 4 . 1 :Figure 4 . 2 :

 4142 Figure 4.1: Feng et al. tested fillers [1] : (a) Metallic foam, (b) Array of fins and (c) Finned foam.

(4. 2 )Figure 4 . 3 :

 243 Figure 4.3: Hybrid filler comparison with existing solutions using: (a) LTE assumption (average and maximum die temperature) or (b) Local Thermal Non-Equilibrium (LTNE) assumption (average and maximum die temperature).

4. 3 . 1 Figure 4 . 4 :

 3144 Figure 4.4: Preliminary study optimization domain.

Figure 4 . 6 :

 46 Figure 4.6: Space colonization algorithm principle.

Figure 4 . 9 :

 49 Figure 4.9: Genetic algorithm flow chart.

CHAPTER 4 .

 4 BEYOND HOMOGENEOUS FILLER generations, a loss of diversity.

Figure 4 . 11 :

 411 Figure 4.11: Reference parameter resulting (a) Network and (b) Heat spreader.

ure 4 .

 4 12e) for the positioning of the auxin nodes.

Figure 4 . 12 :

 412 Figure 4.12: Evolution of the resulting network from the reference one when: (a) D = L X /240, (b) N source = 13, (c) N auxins = 21, (d) r max = L x /60, (e) Half-circles undergo a 180°rotation, (f) N circle = 8,

Figure 4 .

 4 Figure 4.13: Resulting (a) maximum and (b) average fitness for different population sizes over 50 generations running the GDA.

Fitness

  brated. Actually, in part 4.3.2 we discussed the influence of the different parameters influencing the convergence and convergence rate of the genetic algorithm. In this study, we focused mainly on two parameters for time saving reasons: population size N and total generation number M. As a reminder, the population size represents the number of individual initially generated among which best individuals are selected and combined to iteratively improved over M steps.Although, it has been proven that p Off and p Mut also influence convergence and convergence rate, we decided to fix them at respective values of 0.2 and 0.25, which means that 20% of the individuals are replaced at each generation and their replacements see up to 25% of their genes mutated after they have been created through the crossover step. Another parameter was added following Koumousis et al.[19] who proved that introducing random individuals at some generations could improve convergence and convergence rate. At theses specific generations, an additional proportion p Rep of the individuals are eliminated during the selection step and replaced by random ones during the crossover.Those random individuals do not undergo any kind of mutation. It was decided to introduce those random individuals every five generations.We chose to run the GDA for six different population sizes and inspect the evolution of the fitness (maximum and average) over 50 generations (see figure4.13). While the maximum fitness gives us information on convergence and convergence rate, the average one After 5 generations (Simulation data) After 30 generations (Simulation data) After 50 generations (Simulation data) After 5 generations (Curve fitting) After 30 generations (Curve fitting) After 50 generations (Curve fitting)

Figure 4 . 14 :

 414 Figure 4.14: Resulting maximum fitness for different population sizes.

Figure 4 . 15 :

 415 Figure 4.15: Heat spreaders generated by GDA after 50 generations using population sizes of: (a) 5 individuals, (b) 10 individuals, (c) 15 individuals, (d) 20 individuals, (e) 30 individuals and (f) 50 individuals.

Figure 4 . 17 :Figure 4 . 18 :

 417418 Figure 4.17: GDA resulting structure compared temperature distribution and average temperature (T avg ) evolution over time.

Figure 4 . 19 :

 419 Figure 4.19: Use-case optimization space for GDA.

Figure 4 . 20 :

 420 Figure 4.20: Reference heat spreader generated in the use-case optimization domain.

  3.2, source node a regularly spread along the base cylinder limits (see figure 4.20).

Figure 4 . 21 :

 421 Figure 4.21: Evolution of the resulting network from the reference one when: (a) R 0 = 5, (b) N source = 8, (c) Half-circles undergo a 180°rotation, (d) R Source = 2R 0 , (e) N auxins = 20 , (f) N circle = 10, (g) η x = 1.75 and (h) η y = 1.75

ure 4 .

 4 20) for the positioning of the auxin nodes.

Figure 4 . 22 :

 422 Figure 4.22: Spreader growth adaptation to avoid fusion front constraint.

figure 4 .Figure 4 .

 44 figure 4.19 which consists in a 2D projection of the upper plan of the actual LHTESS. As it is done in the final system, the heat source is represented by the sub-domain Γ h while the outer boundaries of the system represented by Γ n are isolated. The change in dimensions

Figure 4 . 26 :

 426 Figure 4.26: Compared top cover influence on the end-of-mission temperature for spreader free or optimized spreader based devices.

3 . 3 ,

 33 finding an optimum individual could require larger population size or generations to cope with the additional parameters necessary to generate the spreaders.

Figure 4 . 27 :

 427 Figure 4.27: Compared average die temperature profiles of devices using hybrid (optimized and non-optimized) and lattice fillers.

- 1 .K - 1 ] 4 Hr

 114 ×10 = 1K/min, m = 3.09 mg Hr = 5K/min, m = 3.09 mg Hr = 10K/min, m = 3.09 mg Hr = 1K/min, m = 10.46 mg Hr = 5K/min, m = 10.46 mg Hr = 10K/min, m = 10.46 mg Isothermal step method Pseudo-Voigt Model fit

Figure A. 3 :

 3 Figure A.3: Influence of the sample mass on the resulting heat flow measurement for two Rubitherm RT80HC samples of respective mass m 1 = 4.02 mg and m 2 = 9.18 mg.

- 1 .K - 1 ] 4 Hr

 114 ×10 = 1K/min, m = 3.22 mg Hr = 5K/min, m = 3.22 mg Hr = 10K/min, m = 3.22 mg Hr = 1K/min, m = 7.21 mg Hr = 5K/min, m = 7.21 mg Hr = 10K/min, m = 7.21 mg Isothermal step method Pseudo-Voigt fit

Figure A. 4 :- 1 .K - 1 ] 4 Hr

 4114 Figure A.4: Influence of the sample mass on the resulting heat flow measurement for two Rubitherm RT90HC samples of respective mass m 1 = 3.83 mg and m 2 = 8.52 mg.

Figure A. 5 :Figure A. 6 :Figure A. 7 :Figure A. 8 :

 5678 Figure A.5: Influence of the sample mass on the resulting heat flow measurement for two Rubitherm RT100HC samples of respective mass m 1 = 3.83 mg and m 2 = 8.52 mg.

  1.1 PCM categories properties. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1 Dimensionless parameters expression and values. . . . . . . . . . . . . . . . 80 2.2 Experimental Validation Sample Dimensions. . . . . . . . . . . . . . . . . . .

2.3 Base Material Measured Specific Heat Capacity. . . . . . . . . . . . . . . . . . 2.4 Base Material Measured Thermal Diffusivity and Conductivity . . . . . . . . 2.5 Samples ETC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Selected PCM melting temperatures, latent heat of fusion and mass of the samples used for the DSC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Thermal properties of the different materials used in simulation. . . . . . . . 3.3 Geometrical parameters of the different PCM based devices. . . . . . . . . . . 3.4 ETC Matrix and specific area of the manufactured PCM based devices . . . . 3.5 Material used in the experimental setup adapted simulation (at 20 • C) . . . 3.6 Experimental and Estimated values of the integral heat transfer coefficient comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 Photos of different samples used to experimentally validate ETC detailed in chapter 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V A.2 Lattice Structures Thermal Diffusivity and Conductivity Measurements and Comparison with both constriction-free and constriction ETC Models. . . . X A.3 Resulting GDA based structures growth parameters. . . . . . . . . . . . . . . . XI

Table 1 .

 1 1: PCM categories properties.

	No	Yes	No	++

  Experimental data using Copper foams(Yao et al.) Experimental data using Aluminum foams(Ranut et al.) Experimental data using Carbon foams(Loeblein et al. and Ji et al.)

	0.1
	0.09
	0.08	Series Model
	0.07	Parallel Model Singh & Kasana
	0.04 0.05 0.06 k e, =k s	Calmidi & Mahajan Yang et al. Boomsma and Poulikakos (Corrected by Dai et al.
	0.03
	0.02
	0.01
	0
	0.85	0.9	0.95	1
		Porosity

2. While correcting Boomsma and Poulikakos, Dai et al. proved that this condition can only be met for high porosity (ε > 0.97) if e is kept constant. Eventually, the issue was discussed by Yang et al.

[8] who pointed out that according to equation 1.79 a decrease in porosity meant an increase in strut diameter (d increases) while e and the node edge remain constant, leading CHAPTER 1. PHASE CHANGE MATERIALS ) Figure 1.24: Foam ETC analytical models comparison. to a geometrical impossibility. To solve the issue, Yang et al. used experimental data to fit e by a third degree polynomial expression, function of the porosity:

5% 93.5% 95.5% 97.5% 88.5% 90.5% 92.5% 94.5% 85.5% 87.5% 89.5% 91.5% 82.5% 84.5% 86.5% 88.5%

  Zhenget al. then observed the evolution of the melting front as well as the thermal energy storage (energy absorbed by the system at a given time) and the thermal energy storage rate.

	Three types of porosity gradients were studied (see figure 1.29): horizontal, vertical and
	two-dimensional; with different levels of gradient and were then compared to a compos-
	ite filled with a uniform foam of similar global porosity. Selecting the optimum porosity
	for each gradient type (see figure 1.29) allows to reduce the total melting time for a given

boundary condition on the left side of the enclosure. From 1320 s for a uniform porosity it goes down to around 1305 s for the optimum horizontal gradient, 1245 s for the optimum vertical gradient and 1185 s for the optimum two dimensional gradient. Before discussing the results, it should be noted that, in this case, heat transfer is mainly convec-91.

  2 × P s in BCC or BCC z cells while L = 2/2 × P s in FCC and FCC z cells. From those CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

				Aluminum
				PCM
	Heat spreaders
	Figure 2.7:
	statements, we draw:
	e =	2 β	πα 3	,	(2.63)
	for a cubic node in a BCC or BCC z cell,
	e =		2πα β	,	(2.64)
	for a cubic node in a FCC or FCC z cell, and
	e =	απ 2β	,		(2.65)

ETC measurement simulation setup.

Table 2 .

 2 3.3. We measure the difference of

	Type	α 1	α 2		e 1		e 2	ω ε	ω z	ω x y
	Foam [22]	3/2	-	see eq. 1.95		-	-	-	-
	BCC	4 π 1.9933 2/3	-	2 β	πα 1 3		-	1.1269 2.3614	-
	FCC	4 π 1.9933 2/3 4 π 0.9041 1/2		2πα 1 β		α 2 π 2β	1.2575 2.3451 2.1892
	BCC z	4 π 3.888 2/3	4 π 1.9933 2/3	2 β	πα 1 3	2 β	πα 2 3	1.0946 2.0033 1.9369
	FCC z	4 π 3.888 2/3	4 π 0.9041 1/2		2πα 1 β		α 2 π 2β	1.4044 1.5224 1.5059

1: Dimensionless parameters expression and values.

  CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING1) The heat flux is specified over S 1 which means that this is not a Mixed Boundary Value Problem (MBVP) (i.e. a problem in which both Neumann and Dirichlet boundary conditions are encountered). This problem is dealt with by Negus

Two

important cases are encountered:

  1 A s /π (i.e. radius of the node in contact with the heat spreader) in a squared compound disk of side P s and height h (see figure2.16). Using results from Muzychka et al.we conclude that it is equivalent as calculating the constriction resistance when heat flow through a contact spot of radius α 1 r in a circular compound disk of radius P

	CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND
	MODELING
	P s /√π
	h	P s
	√α 1 ×r	h	√α 1 ×r
	Figure 2.16: Equivalence of thermal constriction resistance between a circular (left) and a squared
	(right) compound disk.

Going back to the issue presented in section 2.2.3, we need to find a way to calculate the constriction resistance appearing when heat flow through a contact spot of radius α s / π and height h. Substituting in equations 2.102 and 2.103, we can re-calculate θ and φ as:

  over the one defined in equation 2.7 and theoretically validate it, a new series of simulations is run. In that series a single FCC cell surrounded by PCM with 10 µm thick heat spreader on its top and at its bottom is simulated. The simulation is done for four different levels of porosity (i.e. 0.85, 0.91, 0.95 and 0.98) and pore size varying from 200 µm to 3 mm. The fact that ETC is measured

on a unique cell with relatively thin heat spreaders should emphasize thermal constriction making it more visible. We determine the ETC using a similar method to what was done in part 2.2.2 using equation 2.66. We depict the results of the simulations as well as the value given by the different models in figure

2

.

18

. Similar results for the other cell types can be found in Appendix A.1.1. On this figure we notice an almost perfect fit of the model in equation 2.109 with simulations data with deviation for low pore sizes mainly due to numerical limitations. On the other hand, we observe that the first derived model (equation 2.7) is too limited to account for the changes. Accounting only for topology to determine the ETC of such configuration does not seem to be sufficient and accounting

Table 2 .

 2 3: Base Material Measured Specific Heat Capacity.

	CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND
	MODELING
	Temperature ( • C) Specific Heat Capacity (J•kg -1 •K -1 ) Uncertainty (J•kg -1 •K -1 )
	25	884	44
	50	899	45
	100	936	47
	150	968	48
	200	1029	51

.110) Figure 2.23: Enthalpic graph obtained for the base material at 200 • C.

  • C, 50 • C, 100 • C, 150 • C, 200 • C to observe the sensitivity of the specific heat capacity to the temperature. The values and the uncertainty on the measure (5%) are reported in table 2.3. We notice a variation of the specific heat capacity with the temperature which is expected. The values are coherent with those found in the literature for such alloy (C p ≈ 950 J.kg -1 .K -1 at 25 • C).

Table 2 .

 2

		61.1	3.1	137.3
	50	61.3	3.1	140.9
	100	61.2	3.1	145.6
	150	60.2	3.0	148.1
	200	61.1	3.1	159.8

4: Base Material Measured Thermal Diffusivity and Conductivity CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND MODELING

Table 2 .

 2

	CHAPTER 2. PHASE CHANGE MATERIAL ENHANCERS: TOPOLOGY INFLUENCE AND
	MODELING
	Sample	ETC Measured (W.K -1 .m -1 )	ETC Model (W.K -1 .m -1 )	Relative error (%)
			Eq. 2.7 Eq. 2.109 Eq. 2.7 Eq. 2.109
	1	16.86	18.00	17.18	6.31	1.84
	2	16.34	18.00	16.61	9.23	1.63
	3	16.22	18.00	15.92	9.87	1.91
	4	19.66	19.34	18.39	1.48	6.73
	5	37.48	34.58	33.22	8.38	12.81
	6	8.85	8.55	8.03	3.46	10.16
	7	23.56	23.68	22.25	0.49	5.90
	8	22.83	23.68	22.66	3.56	0.78
	9	23.48	23.68	22.78	0.83	3.09
	10	22.58	23.17	22.19	2.53	1.78
	11	22.41	23.17	22.01	3.27	1.84
	12	23.21	23.17	22.33	0.17	3.94

.122)

Looking at the values of the Thermal Diffusivity in table A.2, we notice that this parameter is almost constant on the temperature range considered in this study which is consistent with the results obtained on the reference sample studied in part 2.4.2. 5: Samples ETC.

Table 3 .

 3 described in the previous chapter, such measurement consists in comparing the heat flux inside two different aluminum crucibles while they are exposed to a temperature 1: Selected PCM melting temperatures, latent heat of fusion and mass of the samples used for the DSC.

		Melting	Latent heat	Sample 1#	Sample 2#
	PCM	temperature ( • C)	of fusion (J•g -1 )	mass (mg)	mass (mg)
	Rubitherm RT80HC	78	220	3.09	10.46
	Rubitherm RT90HC	90	170	3.22	7.21
	Rubitherm RT100HC	100	180	4.02	9.18
	Puretemp 108	108	180	3.83	8.52

  • C for a 9 mm thick bare aluminum sample, 52 • C for a 10 mm one and 48 • C for a 11 mm sample, representing in most case around 20% temperature gain. Although this gain is reasonable, it still barely reaches the 120 • C target required for this use-case and different solutions need to be explored. Concept like Cascaded storage, defined in part 1.4.3, should be examined next. CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE Those results, although satisfying, are based on several approximations and assumptions : LTE and absence of convection. As convection is buoyancy induced (see part 1.3.3)

Table 3 .

 3 4: ETC Matrix and specific area of the manufactured PCM based devices

	• Overall height H,
	• Top cover thickness h t ,

Table 3 .

 3 5: Material used in the experimental setup adapted simulation (at 20 • C) ples with respective thicknesses of 9, 10 and 11 mm. Using such sample simplify the simulation as it avoid any melting or filler/PCM relation issue letting us concentrate on the boundary conditions only.In a first attempt to calibrate the simulation, the heat source is no longer represented as a heated surface. Instead, a representation of the actual heating resistance stack-up is represented based on the manufacturer information. As a consequence, it is represented as a three layers stack-up, the first one (on top) being the resisting and heating element made of nickel and chromium which is deposited on an alumina substrate which is itself deposited on a nickel platted copper flange which serves as the interface material with the sample.

	•kg -1 ]

Table 3 .

 3 •K -1 ) 6: Experimental and Estimated values of the integral heat transfer coefficient comparison.

		Measurement	Model	Pore size (mm)	Porosity	Overall thickness (mm)	Top cover thickness (mm)
	1	6.32 × 10 5	4.26 × 10 5	1.75	0.8	9	4.5
	2	3.58 × 10 5	2.58 × 10 5	2.25	0.8	10	4.5
	3	9.05 × 10 5	3.90 × 10 5	1.83	0.8	11	4.5
	4	5.75 × 10 5	2.50 × 10 5	2.75	0.7	11	4.5
	5	5.60 × 10 5	1.73 × 10 5	2.75	0.8	11	4.5
	6	5.40 × 10 5	1.07 × 10 5	2.75	0.9	11	4.5
	7	2.37 × 10 5	8.16 × 10 4	4.00	0.8	11	2
	8	2.00 × 10 5	4.32 × 10 4	5.50	0.8	11	4.5

  .11) Values obtained using experimentally and using this model are reported in table 3.6.This table shows that while values are obviously different as topologies are different, the variations are coherent showing the validity of the experimental values. Actually, while some values are harder to compare, looking at samples #3, #5 and #8, it shows that in both cases h actually depends on the pore size and is roughly proportional to 1/P 2 s .

  3.2 on a device similar to sample #3 with LTE and LTNE assumptions. Results are depicted in figure 3.20 and shows a temperature increase of about 7 • C between LTE and LTNE which places such solution above the 120 • C target temperature. This illustrates how thermal lagging can ac-CHAPTER 3. PCM ARCHITECTED FILLER IMPLEMENTATION ON A USE-CASE tually slow down fusion process, and at the same time energy storing, causing the quick temperature rise. Compared to the 10 • C gain found in part 3.3.2, this limits the interest of PCM based solution for this particular use-case apart from the 24% mass gain.

  the previous chapter, different modeling tools were implemented in a simulation model to assess the performance and optimize aluminum lattice/Phase Change Material (PCM) based Latent Heat Thermal Energy Storage System (LHTESS). Geometry parametric analysis based on these simulations help to determine a number of geometrical parameters that could be significant in enhancing the PCM based device. However, through experiments, it was shown that assumptions made in the simulations tended to overestimate the performances of the devices although the qualitative results remained valid. Anyway, further improvement of both composite and overall system design are required. It was shown, in chapter 1, that further improvement of PCM based systems could be done by combining different types of filler in terms of porosity, topology or even with different PCM in the same system in a process called Cascaded storage.To finally reach the 120 • C target for the use-case defined in part 3.3.1 further PCM enhancement is necessary. As a consequence, in the first part of this chapter and based

on Feng et al.

[1]

, Guo et al.

[2]

Yang et al.

[3]

Finned Foam composite, we develop a new

  .4. The domain Ω is a square of side length L x while the heat source is represented on a Neumann boundary Γ a of length l h . The remaining of the domain boundary Γ n is isolated. Any part of the domain Ω is to be filled either by bare alu-

	CHAPTER 4. BEYOND HOMOGENEOUS FILLER
	Γ h		Heat Spreader
			Enhanced PCM
		[k s , C p,s , ρ s ]
	Ω	[k xy (ε h ), C p,eff (ε h ), ρ eff (ε h )]
	Γ n
	Figure 4.5: Preliminary study domain filled with a heat spreader and enhanced PCM.
	equation:
	k(x)∇ 2 T(x, t ) = C p (x)ρ(x)	∂T ∂t x,t	,	(4.4)

minum that should form a unique shape of area A b or an effective material formed by the combination of PCM and a lattice filler of porosity ε h given by equation 4.2 in which the value of ε 0 determines the quantity of PCM to be found globally in Ω. Once the material distribution is determined, the temperature distribution is evaluated using the following with x ∈ Ω a given point of Ω, the value of thermal conductivity, specific heat capacity and

density are equal to [k s , C p,s , ρ s ] if x is located in the heat spreader [k x y (ε h ), C p,eff (ε h ), ρ eff (ε h )] (

see equations 2.59, 1.36 and 1.35) if it is located in the enhanced PCM (see figure 4.5). For simplification purposes, LTE assumption is considered.

RÉSUMÉ EN FRANÇAISproduite par fabrication additive. Les résultats montrent qu'il est absolument nécessaire de considérer les effet du Thermal Lagging dans les simulations afin d'obtenir une bonne concordance avec les mesures. En utilisant cette hypothèse, nous sommes capables de rendre compte du comportement des dispositifs par des simulations et de valider les résultats de l'étude paramétrique.En partant du constat que la proportion globale de matériau conducteur, mais également son emplacement dans le dispositif, avaient une influence sur la température de la source de chaleur nous proposons, dans le dernier chapitre, d'utiliser la méthode de Cascaded storage. Cette méthode vise à modifier localement la proportion de PCM tout en gardant une proportion globale de PCM constante afin de faciliter la dissipation de chaleur, en maintenant une capacité de stockage équivalente. Pour cela, nous développons un algorithme (Generative Design Algorithm ou GDA) permettant de générer des structures semblables à celle des veines de feuilles particulièrement adaptées à la dissipation thermique. Celles-ci permettent de dissiper efficacement la chaleur proche de la source de chaleur et dans l'ensemble du dispositif homogénéisant ainsi la température, ce qui globalement doit permettre un stockage plus efficace. Le design de la structure est optimisé en utilisant un algorithme génétique. Les résultats des simulations montrent que de telles structures permettent d'améliorer les performances des dispositifs du troisième chapitre.Nous avons prouvé dans cette thèse, que des dispositifs à base de PCM peuvent être une alternative viable aux dissipateurs métalliques classiques. Leur faible densité combinée à leur forte chaleur latente permettent de gagner jusqu'à 20% en masse et en échauffement. Cela a été rendu possible par le développement de nouveaux renforts produits par fabrication additive permettant une meilleure intégration au dispositif et par la réalisation du dispositif suivant la théorie du Cascaded storage. Il est à noter que des progrès peuvent être faits en améliorant les algorithmes qui adaptent la quantité locale de renfort et en améliorant les procédés de fabrication additive, ce qui permettra de réaliser des structures à plus petites échelles et de réduire les effets de Thermal Lagging et de constriction. viii

see Rubitherm GmbH or PureTemp for examples.

Acknowledgments

ACKNOWLEDGMENTS

x

Base Material Characterization

Before the samples depicted in Appendix A.1.2 can be studied, a precise knowledge of the material they are made of is necessary. Actually, Flash Laser method only provides information on the thermal diffusivity of the studied sample which means that a precise knowledge of both specific heat capacity (C p ) and density (ρ) is required. To compare the experimental data to the model, a precise determination of the thermal conductivity of the base material is also necessary as it is used to calculate the ETC. Due to the complexity of the measurements combined with necessity of using precise machines that

were not available at IEMN, it was decided that the measurement would be performed by Influtherm, a company specialized in thermal related measurement.

Base Material properties measurement

As mentioned in the previous paragraph a precise knowledge of a number of properties is necessary to derive the thermal conductivity from the thermal diffusivity. In that sense, three intrinsic properties of the samples base material are measured : density (ρ), specific heat capacity (C p ) and coefficient of thermal expansion (α) that may have an impact on the diffusivity measurement.

The first property to be measured is the density. For that purpose a cuboid made of the base material is used which dimensions, measured using a caliper, are as follows:

• length: 11. 64 • Pore size P s ,

• Number of stacked cells N,

• Porosity ε.

In addition to the parameters already described in the previous paragraph two additional parameters were added: Pore size and Number of stacked cells. Although these parameters were neglected in the previous paragraph, it was already shown in chapter 2 that both could affect the ETC. In addition, the pore size also affects the specific area that is related to the LTE assumption. To illustrate both facts, we report in table 3.4 the value of the ETC matrix and specific area a v . The specific area is calculated using the following expression using the parameters defined in part 2.2.2:

(3.9)

We notice a slight variability of the ETC in the vertical direction for similar porosity due to thermal constriction which has been linked to both pore size as well as bottom and top cover thicknesses (see part 2.3.2). On the other hand, the specific area depends on both pore size and porosity as seen in equation 3.9.

These specific designs allow us to assess the influence of :

• combined pore size and number of stacked cells comparing samples #3, #5 and #8,

• combined pore size and overall package thickness comparing samples #1, #2 and #5,

• porosity comparing samples #4, #5 and #6,

• combined pore size and top cover thickness comparing samples #5 and #7.

Chapter 4

Beyond Homogeneous Filler

The art of structure is where and how to put the holes. Using the iterative process, N growth new vein nodes are created and their position is determined by:

Robert le Ricolais

with D the spatial step between each vein node and k represents the set of all growth directions at the step i and n = n 1 . . .

. The algorithm is run recursively until all auxin nodes have been reached (i.e. a vein node is found at a distance inferior to the kill distance d k ). An example of such process is found in figure 4.6.

Once the network growth is completed, each branch is selected and sorted from longest to shortest. Then, extremity node radius is set to a minimum value r min and radius grows linearly until the original source is reached and which radius is set to the maximum value r max . The longest branch of length L b serves as a reference to define linear coefficient of radius variation:

From it, each node radius is defined as a function of its distance to the branch extremity d :

Another algorithm is developed to determined the boundaries of the resulting shape.

From this new set of points and the use of an alpha shape allow the definition of the final shape. An alpha shape is one representation of a shape that contains a given set of points.

Actually, a "shape" has no formal meaning and it is especially true when we are looking for a "shape containing a given set of points". To define an alpha shape, we consider a given domain of R 2 and a given set of points σ. For each point of σ we define, if possible a circle of radius r α that does not enclose any of the other point of σ. The alpha shape is given by the set of points S p at which it was possible to create such circle (see figure 4.8) [14]. A variety of shapes can be obtained for a given set of points just by changing the value of r α .

In our case we are setting it as the minimum value that gives a single region for our set of points leading to the creation of the shape in figure 4.7. One of the main advantage of this generation process is that the final shape, rather than being describe entirely, can only be described by a set of parameters determining the position of the source and the auxin nodes and the minimum and maximum network node radius. Using such representation, the set of parameters can be optimized using an simple and robust algorithm, particularly adapted to low dimension optimization: the genetic algorithm due to Holland [15].

Genetic algorithm

The space colonization algorithm is used to generate heat spreaders using a limited number of parameters that entirely describe its shape. This way of describing the heat spreader greatly reduces the dimension of the design problem and will represent the vector s defined in equation 4.3. Using a biological analogy, the vector s including the different parameters describing the structure is called the genotype and the resulting heat spreader shape the phenotype. Following the analogy, the theory of Charles Darwin on natural evolution stating that natural selection occurs for "the fittest to survive" is used to select the best structure (or individual) and use them as a way to produce new structures (or offsprings) inheriting some parts of their phenotype. Over a certain number of generations (or iteration of the process), we expect the maximum and the average fitness of the population to increase.

Practically , the vector s described earlier contains a series of parameters that can be of four types [17]:

• Binary,

• Real number,

• Integer/literal permutation,

• A combination of the three previous types. Each of the parameter is called a chromosome that can take different values (of a given type for each chromosome) called alleles [18]. The use of one type or the other is problem dependent but in many real-life problem a combination of different types is often necessary.

At the initialization step of the algorithm a series of N different individuals, each represented by a unique vector s i , is created forming a so-called population. The allele of each of their chromosome is randomly selected in a attempt to cover as many part of the optimization space (representing the whole of the possible individuals) as possible. During the second step of the algorithm, all individuals of the population are tested and their fitness function F is evaluated. The way of testing each individual highly depends on the problem, though in our case it is done through FEA. We are sensing here that while larger values of N provides a higher genetic diversity and a better exploration of the optimization space, it also leads to higher calculation times which is detrimental to the algorithm performances: in that sense an optimum value for N is to be found for each problem depending mainly on the number of chromosomes and the time necessary for the evaluation.

Once all individuals of the population have been evaluated, the algorithm selects the Comparing the different hybrid fillers, we notice that the five worst individuals behave in a similar way in terms of both heat source temperature and temperature increase rate whereas small differences are noticed in terms of temperature distribution and average temperature (see figure 4.17). The best individual (obtained with a population size of 50)

behave, on the other hand quite differently. Heat source temperature increases less at the start of the simulation. This fact could be explained by the smaller surface A b occupied by the sixth spreader: while every other spreaders occupy more than 15% of the total domain surface this one only occupies around 10% of it, leaving more conductive material to be used as filler and increasing the thermal conductivity of the enhanced PCM. This provide an easier path for the heat to spread in the entire structure: close to the heat source, the high diffusivity material helps avoiding any energy concentration while heat can still efficiently spread in the enhanced PCM where it is stored. As seen in figure 4.17, in this case the temperature is more homogeneous in the entire structure. This does not necessarily increase the melting rate (as seen on figure 4.18), which is usually associated with better heat storage in the literature, but the entirety of the LHTESS can be better used during the entire simulation duration. A better distribution of conductive material inside the entire system, locally reduces thermal gradient, avoiding hot spot that, in this particular case, is always located at the heat source.

Looking more specifically at the temperature profile of the heat source using the best obviously affects the behavior of the system: while in the 3D system, the heat source is represented as a heat flux entering a boundary of the system, in a domain made of aluminum, in this case, the heat source has the same dimensions as the system itself and is applied in sub-domains made of both aluminum and enhanced PCM which, following Fourier's law might induced a higher temperature gradient in the sub-domains made of the latter. In addition, heat transfer is not facilitated by the top cover. For these reasons, the value of power input are reduced to 20 W instead of 100 W and the time range increased to 700 s instead of 30 s. Although, this obviously affects the results of the simulation, it should give us insights on whether a given structure is able to efficiently spread heat in the x y plan to maximize melting rate and heat storage or not.

To determine an optimal spreader, the GDA is run with three different initial population of fifty individuals over one hundred generations with similar offspring and mutation probabilities than in part 4.3.3. The resulting maximum and average fitness for each generation and population are reported respectively in figures 4.23a and 4.23b. The latter shows that the genetic algorithm is able to maintain a genetic diversity over the generations as the average fitness tends to vary rather than constantly increase. On the other hand, looking at maximum fitness over generations, it seems that a close-to-optimum solution is found in run number 3 as no progress is made over more than 70 generations.

In figure 4.24, we report both temperature and temperature increase rate profiles for the three different optimized spreaders. Comparing them to the fusion front depicted in figure 4.25, it appears once again that any constraint to the fusion front propagation leads In addition, a faster melting rate is reached for spreader based devices as total melting occurs sooner.

Transition to 3D and discussion

The previous 2D study helped determining three optimized heat spreader to create a hybrid filler that might solve the issues discussed in chapter 3 related to LTNE and more globally efficient heat spreading. Although the GDA gave decent results allowing a 5 • C gain compared to homogeneous filler, it was only done in 2D for time constraints. This completely modifies boundary conditions which may lessen the efficiency of the optimized spreaders.

To confirm or infirm this last assumption and determine the actual ability of the optimized structures to create an efficient hybrid filler, they are integrated in an actual usecase, like it was done for simplified heat spreader in part 4. profiles of these devices with the one obtained in chapter 3 using a homogeneous filler, we showed that the hybrid filler device outperformed the two others which was attributed to the combined effect of efficient heat spreading close to the heat source and decent melting rate insured by the enhanced PCM.

From the results obtained with simplified heat spreader, it was then decided to find a way to optimize the spreader topology in order to further increase its performances.

To do so, a so-called generative design algorithm (GDA) was developed. This algorithm combines a space-colonization algorithm for the spreader generation given a set of parameters and a genetic algorithm for the optimization of those parameters. To assess the ability of this algorithm to enhance heat transfer and determine the suitable parameters, it is first used on a 2D case and the resulting structures are compared to an homogeneous filler showing encouraging results. This preliminary study also helped determine the suitable population size and number of generations in order to maximize fitness and convergence rate. It was concluded that while the number of generations could be limited, increasing the population size was best to prevent premature convergence and while improving convergence rate.

In the final part of the chapter, the GDA is used on the actual use-case, first on a simplified 2D version to save simulation time and optimize the structure heat transfer performance in the x y plan. The 2D optimization shows an up-to 5 • C gain compared with homogeneous filler solution. The three resulting optimized spreaders are then integrated into the actual device to assess their respective performances under both LTE and LTNE assumptions. The simulations shows very similar results for the three optimized spreaders using both assumptions which was expected due to their behavior in 2D. Finally, those spreaders, optimized in 2D, are compared with the hybrid filler defined in partı4.2. It showed that the four spreader behave in a similar way, showing that additional optimization strategies need to be employed.

[1] S. Feng, M. Shi, Y. Li, and T. J. Lu, "Pore-scale and volume-averaged numerical simulations of melting phase change heat transfer in finned metal foam," )

Temperature ( Phase Change Material (PCM) have been widely used for thermal energy storage due to their high latent heat of fusion for thermal management applications. However, PCM suffer from their very low thermal conductivity which limits heat spreading around the heat source. Without proper thermal conductivity enhancement, melting would occur mainly around the heat source and heat would be conducted too slowly for the device to be efficient. To do so, a variety of fillers have been Les matériaux à changement de phase (PCM) sont largement utilisés du fait de leur forte chaleur latente de fusion pour des applications de management thermique. Malgré cela, les PCM pâtissent de leur faible conductivité thermique qui limite la propagation de la chaleur autour de la source.

Sans amélioration notable de la conductivité, la fusion ne se déroule qu'autour de la source de chaleur limitant grandement les performances du dispositif. Pour parer à cela, une multitudes de renforts thermiques ont été développées et étudiées ces dernieres années. Une revue de littérature détaille l'étude de ces différents renforts ainsi que de leurs performances lorsqu'ils sont combinés aux PCM dans le premier chapitre de ce manuscrit. Le deuxième chapitre présente, lui, un nouveau renfort basé sur l'utilisation de lattices métalliques ainsi qu'un modèle permettant