SIMULATION NUMÉRIQUE ET QUANTIFICATION D'INCERTITUDES POUR LE REFROIDISSEMENT

Keywords: Batteries Lithium-ion, Refroidissement par immersion, Simulation numérique, Quantification d'incertitudes, Transfert thermique conjugué, Calibration Bayésienne, Métamodèles Lithium-ion batteries, Immersion cooling, Numerical simulation, Conjugate heat transfer, Uncertainty quantification, Bayesian inverse problem, Surrogate models

PAR IMMERSION DES BATTERIES LITHIUM-ION

Résumé Pour promouvoir une utilisation plus large des véhicules électriques, les batteries Lithium-ion (Li-ion) se doivent de supporter des courants électriques importants, générant ainsi de fortes contraintes thermiques qui dégradent leurs performances et leur durée de vie. La gestion thermique des packs de batteries est donc un élément crucial pour répondre à ces nouvelles contraintes industrielles. La technologie de refroidissement par immersion est une solution prometteuse en termes de performances thermiques. Ces systèmes sont gouvernés par des phénomènes multi-physiques, allant de la chimie interne des batteries jusqu'au transfert thermique à l'échelle du pack de batteries alimentant le moteur électrique.

Cette thèse a pour objectif de développer des modèles numériques pour le refroidissement par immersion des batteries Li-ion, tout en considérant les incertitudes provenant des paramètres physiques en jeu.

Cette problématique est abordée en proposant deux modèles de fidélités croissantes. D'abord, un modèle dit basse fidélité est développé, incluant la modélisation des phénomènes thermiques et électriques du problème. Des méthodes de quantification d'incertitudes (calibration Bayésienne et analyse de sensibilité) couplées avec des données expérimentales originales offrent ainsi des éléments de compréhension et d'analyse sur le comportement global du système. Ensuite, une approche plus spécifique est présentée à l'aide d'un code de calcul CFD haute fidélité. Le calcul du transfert thermique conjugué sous régime transitoire en deux dimensions d'un pack de batteries immergées est ainsi réalisé. Cet outil est utilisé pour évaluer la précision d'un modèle construit a priori, représentant la résistance interne des batteries Li-ion. Les incertitudes provenant de la résistance sont prises en compte grâce à la paramétrisation de ce modèle et calibrées en utilisant un cas test expérimental de la littérature. Enfin, pour obtenir une meilleure compréhension de la physique de ces problèmes de refroidissement par immersion, la fidélité de l'outil CFD est augmentée en considérant des calculs de transferts thermiques en 3D, ainsi qu'un modèle de résistance interne amélioré.

REMERCIEMENTS

Il est l'heure d'adresser mes sincères remerciements aux personnes qui ont contribué de près ou de loin à la réussite de cette thèse. De nombreuses personnes ont été impliquées dans les travaux qui vont suivre, sur des plans divers et à des strates différentes. J'ai avant tout vécu ce doctorat comme un enrichissement professionnel, intellectuel et personnel grâce à toutes les personnes que voici. Tout d'abord, mes premiers remerciements se dirigent vers les membres du jury, messieurs Christophe CORRE, Sofiane KHELLADI, Guy FRIEDRICH et Pierre LUBIN qui ont accepté de prendre du temps pour évaluer mes travaux. Je remercie en particulier M. CORRE et M. KHELLADI pour leurs rapports et leurs remarques qui ont permis d'améliorer la qualité de ce manuscrit.

Naturellement, mon deuxième mot de remerciements s'adresse à mes directeurs de thèse, Héloïse BEAUGENDRE et Pietro CONGEDO, pour votre encadrement au cours de ces années.

Merci pour votre aide et vos conseils si précieux, votre disponibilité et pour avoir toujours rendu les choses simples quand elles me paraissaient compliquées. J'ai énormément appris à vos côtés.

Je veux ensuite remercier l'entreprise Exoes, qui sont aussi à l'origine de ce projet de thèse, pour le partenariat très prolifique que nous avons eu pendant ces trois années. En particulier Rémi DACCORD et Maxime GUADAGNINI, pour toute votre expertise sur les batteries et les méthodes de refroidissement. Les réalisations de cette thèse sont issus de travaux que nous avons mené ensemble, avec toujours un vrai intérêt de votre part et des échanges riches de connaissances.

Également, je dois remercier toutes les personnes liées à l'équipe de TrioCFD du CEA Saclay. En particulier Ulrich BIEDER, pour ton aide sur la partie calculs CFD, et pour toutes les discussions passionnantes que nous avons pu avoir au cours de nos nombreuses réunions. J'ai énormément appris à tes côtés et n'aurait pu espérer d'expérience plus enrichissante. Aussi, je tiens à remercier Elie SAIKALI, pour ta contribution sur la partie conductivité anisotrope, qui a permis de commencer à explorer de nouvelles pistes. Enfin, Nabil DJATI pour ton aide et tes conseils techniques sur l'installation et l'utilisation de TrioCFD. J'exprime aussi un grand merci à toutes les personnes du support technique Inria, en particulier l'équipe PlaFRIM, pour leur immense aide sur l'utilisation du cluster.

Je tiens à remercier aussi toutes les personnes de l'équipe CARDAMOM, pour les conseils de chacun et diverses discussions plus générales, et leurs avis extérieurs sur mon travail, toujours pertinents et utiles. Un mot particulier pour notre chef d'équipe, Mario RICCHIUTTO, qui a accepté aussi que je réalise mon stage dans l'équipe et m'a permis de débuter cette aventure à Inria. Également, un grand merci à Anne-Laure GAUTIER pour ton aide lors des multiples démarches et pour toujours avoir su répondre à mes interrogations diverses et variées. D'un point de vue plus personnel je vais remercier tous les PhD-mates, Post-docs et ingénieurs qui iii sont actuellement, ou étaient, dans l'open-space B437, pour toutes les discussions sérieuses et non sérieuses, leur aide et conseils, et surtout pour tous les bons moments passés ensemble, à Inria autour d'un café et dans Bordeaux, autour d'une bière ou de pizzas. Merci à tous.tes et aussi grazie a tutti i chicci pour le petit 'séjour Erasmus' en Italie, it was really ununfun ! Enfin mon dernier mot s'adresse à mes amies et amis de Bordeaux, de Feed The Monkey, d'un peu partout en France et bien sûr, à ma famille. Même si vous ne vous en rendez pas forcément compte, votre présence dans tous ces moments hors du travail a compté beaucoup et m'a aidé grandement à garder la tête froide et reposée pendant ces années. Merci encore, tout ce que j'ai pu réaliser est aussi lié à votre soutien inconditionnel et constant.

RÉSUMÉ ÉTENDU EN FRANÇAIS

Contexte

L'augmentation des émissions de gaz à effets de serre est une menace réelle pour le bien être de l'Humanité sur Terre dans le futur. Dans ce contexte, le secteur des transports est responsable de 30% des émissions totales [START_REF]Ministry for an ecological and solidary transition[END_REF]. Les véhicules électriques propulsés par des batteries Lithium-ion (Li-ion) sont une solution prometteuse pour aller vers des modes de transports à faibles émissions. Les industriels produisent des efforts conséquents pour développer de tels véhicules.

Pour atteindre des parts de marché suffisantes, l'enjeu est de développer des véhicules avec une puissance moteur suffisante et une bonne autonomie. En somme, les véhicules électriques doivent montrer qu'ils sont capables de reproduire les performances des véhicules à propulsion thermique. Ces performances impliquent que les batteries doivent fournir une puissance suffisante sur des longues périodes de fonctionnement. De plus, le temps nécessaire pour une recharge du pack batteries ne doit pas excéder une dizaine de minutes. Les batteries Li-ion sont une solution prometteuse pour atteindre de telles performances, étant donné que ces dispositifs présentent une forte densité de puissance et d'énergie (voir Fig.

II.1)

. Dans les véhicules électriques, des milliers de cellules sont assemblées dans des modules, formant ainsi le pack batterie. C'est ce pack qui fournit l'énergie nécessaire au moteur électrique et fait ainsi avancer le véhicule.

Les performances citées ci-dessus se traduisent en pratique par des courants électriques forts qui sont appliqués aux batteries, pour les séquences de charge ou de décharge. Les réactions électrochimiques produisant l'énergie délivrées par les batteries produisent de la chaleur [START_REF] Sato | Thermal behavior analysis of lithium-ion batteries for electric and hybrid vehicles[END_REF]. Cette génération de chaleur est d'autant plus forte que les courants électriques appliqués sont élevés. Les fortes températures atteintes par les batteries endommagent les batteries à court ou moyen termes, réduisant ainsi leurs durées de vie et leurs performances.

La gestion thermique des batteries pendant leur fonctionnement est donc un élément clé pour atteindre les performances permettant de réaliser des trajets longs à vitesse satisfaisante, toute en réduisant les temps de recharge.

Pour résoudre ces problématiques thermiques liées au fonctionnement des batteries Liion, les industriels tentent de développer des systèmes capables de refroidir les batteries pendant leur fonctionnement. Ce sont les systèmes de gestion thermique des batteries, ou Battery Thermal Management System (BTMS) en anglais. La technique d'intérêt dans le cadre de cette thèse est le refroidissement par immersion. Le rôle de ces systèmes est de (i) maintenir les batteries dans une plage de température garantissant leur sécurité pendant le fonctionnement, (ii) assurer une bonne homogénéité de températures entre toutes les batteries du pack, (iii) assurer une bonne homogénéité de température pour une seule batterie. Dans les systèmes de refroidissement par immersion, les batteries stockées dans un module sont v immergées au contact d'un fluide (de l'air ou un liquide diélectrique). Le fluide circule à travers un circuit de refroidissement, dissipant ainsi la chaleur produite par les batteries pendant leur fonctionnement par transfert thermique convectif.

Simulation des systèmes de refroidissement par immersion pour batteries Li-ion

La simulation des systèmes de refroidissement par immersion, depuis l'échelle du comportement électrochimique interne dans les batteries, jusqu'au transfert thermique entre les batteries et le fluide, consiste en pratique à résoudre un problème multi-physique. Plusieurs approches existent dans la littérature pour réaliser la simulation des BTMS.

Les auteurs dans [START_REF] Shahjalal | A review of thermal management for Li-ion batteries: Prospects, challenges, and issues[END_REF] identifient les méthodes de modélisation permettant de représenter les comportements thermiques et électriques des batteries Li-ion. L'approche usuelle consiste à représenter les batteries et le système de refroidissement avec une géométrie en 1D, 2D ou 3D. Plusieurs types de modèles sont disponibles dans la littérature pour représenter les comportement thermiques et électriques des batteries.

Des modèles 0D et 1D développés dans certaines études ont montré des bonnes capacités à reproduire les prédictions de températures vis à vis de données expérimentales [START_REF] Jaguemont | 1D-Thermal Analysis and Electro-Thermal Modeling of Prismatic-Shape LTO and NMC Batteries[END_REF][START_REF] Amini | A thermal model for Li-ion batteries operating under dynamic conditions[END_REF].

Cependant, l'approche la plus répandue dans la littérature est d'utiliser les outils de la CFD (Computational Fluid Dynamics) pour la simulation des systèmes de refroidissement par immersion. L'approche CFD permet d'obtenir plus d'informations sur le comportement des systèmes, puisque la résolution du système d'équations aux dérivées partielles en 2D ou 3D permet d'obtenir les champs de températures dans le fluide et le solide, ainsi que les champs de vitesse et pression pour le domaine fluide. En se plaçant à l'échelle du pack batterie, le transfert de chaleur entre les batteries et le fluide de refroidissement est usuellement calculé en utilisant l'approche du transfert thermique conjugué. Dans la plupart des études proposant cette approche, les équations de Navier-Stokes sont résolues dans le fluide. L'équation de la chaleur est résolue dans le domaine solide représentant les batteries [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF][START_REF] Lin | A review on recent progress, challenges and perspective of battery thermal management system[END_REF]. Le couplage est réalisé en appliquant la continuité des flux et températures à l'interface.

Les principaux objectifs des études réalisant la simulation des systèmes de refroidissement par immersion peuvent se résumer ainsi: calculer l'évolution de température à des positions spécifiques du pack batteries, suivant des conditions de charge/décharge, estimer le coefficient de transfert thermique pour comparer les performances de refroidissement selon la configuration considérée, le fluide utilisée ou la géométrie de la disposition des cellules dans un pack.

Les études CFD de la littérature sont variées et se focalisent sur différents aspects de la physique des systèmes de refroidissement par immersion. Les études numériques de ces systèmes remplissent un des deux objectifs suivants. Le premier objectif consiste à utiliser les modèles CFD pour caractériser et analyser des processus typiques de transfert thermique dans des cas tests présentant une configuration simplifiée. Le second objectif consiste à utiliser cet outil pour évaluer les performances d'une méthode de refroidissement innovante ou d'un dispositif existant afin de le dimensionner. vi La simulation numérique de ces systèmes présente certaines problématiques. Premièrement, comme vu ci-dessus, les packs de batteries sont composés d'un très grand nombre de cellules, des centaines voire des milliers. Peu d'études réalisent la simulation de si grandes configurations. Représenter les transferts thermiques pour un pack complet peut s'avérer très complexe. De plus, simuler ces configurations de façon exhaustive, en incluant les phénomènes se produisant dans la structure interne des batteries peut s'avérer prohibitif avec des modèles CFD.

Le maillage utilisé pour ces simulations peut être également difficile à réaliser. En effet, les industriels visent à disposer autant de batteries que possible dans un volume donné, rendant l'espace entre les batteries très étroit, parfois inférieur au millimètre. La résolution des équations régissant l'écoulement du fluide peut s'avérer donc très complexe dans ces zones.

De plus, le fluide de refroidissement est en général un fluide complexe, au comportement thermique non classique. Enfin, la chaleur produite par les batteries résulte des réactions chimiques dans la structure interne des batteries. C'est pourquoi, représenter la physique au niveau chimique jusqu'à l'échelle des transferts thermiques au niveau du module demande de développer un modèle multi-physique couplant l'information provenant de ces différentes échelles. Seulement quelques études, par exemple [START_REF] Suresh Patil | A novel dielectric fluid immersion cooling technology for Li-ion battery thermal management[END_REF][START_REF] Jithin | Numerical analysis of single-phase liquid immersion cooling for lithium-ion battery thermal management using different dielectric fluids[END_REF], réalisent de telles approches.

Enfin, les problématiques décrites ci-dessus peuvent expliquer la difficulté à trouver des études avec une description complète du problème multi-physique. La reproduction des résultats de la littérature peut s'avérer parfois difficile en raison du manque d'information disponibles dans ces études sur les divers paramètres physiques mis en jeu.

Quantification d'incertitudes pour le refroidissement par immersion des batteries Li-ion

Les systèmes de refroidissement par immersion présentent un nombre considérable de paramètres incertains à considérer. Premièrement, selon l'échelle considérée, des lois physiques doivent être mises en places pour des paramètres essentiels, régissant les comportements thermique et électrique des batteries.

Plus précisément, les auteurs dans [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF] identifient les principales échelles à considérer pour représenter ces systèmes, en allant des batteries jusqu'au système complet. Premièrement, le comportement interne des batteries peut être résolu au niveau électrochimique.

Il peut être aussi représenté de manière plus macroscopique en considérant les paramètres électriques des batteries. Ici la génération de chaleur est calculée avec des lois phénoménologiques mettant en jeu ces paramètre électriques, comme par exemple l'effet Joule [START_REF] Bernardi | A General Energy Balance for Battery Systems[END_REF]. Enfin, le transfert thermique entre les batteries et le fluide peut être résolu pour une batterie ou au niveau d'un module complet, impliquant de considérer plusieurs batteries et une résolution d'écoulements complexes. Selon l'échelle sélectionnée, des hypothèses doivent être réalisées pour représenter un phénomène physique qui ne sera pas calculé directement.

Toutes ces approches utilisent des modèles plus ou moins distant de la physique interne des batteries, et mettent en jeu des paramètres présentant des incertitudes [START_REF] Der | Aleatory or epistemic? Does it matter?[END_REF]. Autrement dit, il est difficile de choisir des valeurs pour ces paramètres avec un bon niveau de confiance vii si aucune information supplémentaire n'est apportée. Lorsqu'il s'agit d'évaluer les performances des systèmes de refroidissement, les inconnues sur les paramètres vont potentiellement impacter fortement les prédictions obtenues pour les quantités d'intérêts. D'un point de vue de la simulation, ces incertitudes vont impacter fortement les valeurs calculées et prédites par les modèles. La prédiction de ces quantités se montrera plus robuste et informative si la variabilité de ces paramètres d'entrée peut être quantifiée.

Une approche Quantification d'incertitudes permet de résoudre ces problématiques dans le cadre de la simulation numérique. Dans ce cadre, les résultats des simulations ne sont plus vus comme ayant une seule valeur déterministe, mais comme une prédiction avec une probabilité associée qu'il s'agira alors de quantifier. La précision des modèles est alors évaluée en termes de distributions et statistiques (moyenne, variance, quantiles,...) sur les quantités calculées.

Quelques références ont appliqué une approche quantification d'incertitudes pour la simulation des batteries Li-ion. La plupart des études se focalisent sur les équations régissant le comportement électrochimique dans une batterie [START_REF] Vahid Esfahanian | Capacity Fade Analysis of Lithium-Ion Cells Utilizing Uncertainty Quantification Approach[END_REF][START_REF] Hadigol | On Uncertainty Quantification of Lithium-ion Batteries: Application to an LiC 6 LiCoO 2 cell[END_REF][START_REF] Laue | Model-Based Uncertainty Quantification for the Product properties of Lithium-Ion Batteries[END_REF][START_REF] Andrew | Stochastic capacity loss and remaining useful life models for lithium-ion batteries in plug-in hybrid electric vehicles[END_REF]. Quelques auteurs seulement étudient le couplage entre les équations au niveau électronique avec les équations thermiques pour calculer la chaleur produite par les batteries, en considérant les incertitudes à ces deux échelles [START_REF] Tong | Correlating uncertainties of a lithium-ion battery -A Monte Carlo simulation: A Monte Carlo simulation of a lithium-ion battery model[END_REF][START_REF] Chen | Modeling and characterization of the mass transfer and thermal mechanics of the power lithium manganate battery under charging process[END_REF]. Peu de références étudient les effets des incertitudes sur le transfert thermique entre les batteries et le fluide de refroidissement.

Motivations et contributions de la thèse Collaboration industrielle avec Exoes

Dans ce contexte industriel, le projet de cette thèse est issu de la collaboration avec l'entreprise Exoes située à Gradignan en France [START_REF] Exoes | Avenue de la Grande Lande[END_REF]. Exoes développe des modules batteries avec un système de refroidissement par immersion à destination des constructeurs de véhicules électriques. La solution proposée par Exoes est un pack batterie répondant aux besoins des clients en terme d'espace disponible dans le véhicule, de puissance requise par le moteur et de limitations des coûts. Leur expertise est portée sur la construction et l'analyse de ces systèmes de refroidissement.

La collaboration avec l'Inria vient naturellement du besoin exprimé par Exoes de développer des outils numériques avancés pour réaliser des analyses plus poussées de ces systèmes. La complexité des systèmes de refroidissement par immersion nécessite le développement de solveurs haute-fidélité ainsi que d'un traitement rigoureux des différentes source d'incertitudes. De plus, les données expérimentales produites par Exoes permettront une validation systématique des différents modèles numériques.

Motivations

Considérant tous les aspects revus ci-avant, partant de la conception des systèmes de refroidissement par immersion jusqu'au traitement des incertitudes, les objectifs de cette thèse viii peuvent être déclinés ainsi: Simulation numérique des systèmes de refroidissement par immersion pour les batteries Li-ion Le premier objectif de la thèse est de développer des modèles numériques pour réaliser des simulations des systèmes de refroidissement par immersion des batteries Li-ion. Les modèles devront permettre de prédire l'évolution de quantités thermiques comme l'évolution de la température sous différentes conditions. Un compromis entre coûts de calculs et précision sera exploré: deux modèles à différents de niveaux de fidélité seront considérés. La nature des résultats produits par les deux modèles sera mis en perspective avec les ressources de calcul requises.

Identifier et caractériser les incertitudes provenant des systèmes de refroidissement par immersion pour les batteries Li-ion

Ces systèmes multi-physiques nécessitent de traiter les incertitudes afin de réaliser des simulations informatives et prédictives. Cet objectif consiste à produire une analyse physique du processus de refroidissement par immersion pour identifier quels paramètres seront considérés incertains.

Basée sur une revue de la littérature et d'expertises d'ingénieurs d'Exoes, ces paramètres seront sélectionnés. La collaboration avec Exoes permet d'apporter une quantité d'information significative pour caractériser avec précision les incertitudes de ces systèmes. Les incertitudes considérées (plage de variation et distribution) seront définies à ce stade pour les paramètres concernés.

Validation robuste et calibration des modèles vis à vis de données expérimentales Une

fois que les modèle numériques seront mis en place et les incertitudes proprement caractérisées, l'objectif final est de produire des résultats de simulation pour obtenir des éléments d'analyse sur les systèmes de refroidissement par immersion. Tout d'abord, les modèles devront être validés vis à vis de données expérimentales.

L'utilisation directe de ces données permettra ensuite de calibrer les modèle dans un cadre probabiliste (calibration Bayésienne). Les résultats de telles analyses conduiront à des valeurs sur les paramètres incertains, permettant d'obtenir des réponses du modèle proches des données expérimentales. En ce sens, les incertitudes sur les prédictions numériques seront aussi réduites.

Contributions

En suivant les objectifs décrits, les contributions de ces travaux de thèse peuvent se décliner selon les quatre points suivants.

Solveur basse-fidélité à faible coût de calcul

LIST OF ACRONYMS

INTRODUCTION

Chapter abstract

This chapter gives the essential elements to the whole work presented in this thesis manuscript. The industrial and scientific contexts are described to draw the main objectives addressed in this thesis. The main contributions outcoming of this work are summarized. Finally, the content of each chapter composing this manuscript is briefly described so the reader can better understand how the contributions are presented all along the manuscript.

Clean mobility and challenges with Li-ion batteries for Electric Vehicles

The rise of greenhouse gas emissions is a concerning threat to humankind's well-being on Earth in the future. In this context, for instance, in France, the transport activities are responsible for 30% of the total greenhouse gas emissions [START_REF]Ministry for an ecological and solidary transition[END_REF]. More precisely, Fig. I.1 shows that in the EU, the primary source of greenhouse gas emissions is coming from road transport. Electric vehicles (EV) powered by Lithium-ion (Li-ion) batteries stand as a promising solution for zero-emission transportation. Industrials put a lot of efforts on developing such cars: in 2020, the electrical vehicles represented a share of 24.5% of the new cars in EU [START_REF]Making the transition to zero-emission mobility, 2021 progress report. European Automobile Manufacturers' Association[END_REF].

The challenge towards a mass-market adoption of electric vehicles is to develop vehicles with sufficient power and a good range of autonomy. In other words, electric cars should present performances similar to the current internal combustion engine cars. Significantly, the last barrier to large usage of EVs is to prove their capacity of performing long travels across a country [START_REF] Didier | E-mobility, past, present, future. fr. Speech at the NAFEMS France Seminar "Simulation for electric mobility" supported by SIA[END_REF]. It means that the batteries should produce enough power over long periods to allow individuals to travel distances similar to cars with internal combustion engines. As stated in [START_REF] Andre | Future generations of cathode materials: An automotive industry perspective (Review)[END_REF], the coming generations of electrified vehicles would be required to travel 500 kilometers to reach a significant share of the automotive market. Also, the full charging time of the battery pack should not exceed around 15 minutes if one wants to compete with the time required to refill a fuel tank.

Li-ion batteries are the most promising solution to reach such performances, as these devices show excellent energy and power density, being the amount of energy (resp. power) produced by the battery divided by its mass (see Fig. II.1 from chapter II for more details).

Thousands of Li-ion batteries are assembled in modules in electric vehicles, composing the so-called battery pack. This battery pack powers the electric engine, moving the vehicle.

The capabilities required by electric cars stated just above are traduced in practice by high electric currents submitted to the Li-ion batteries for the discharging or charging sequences.

The internal electrochemical reactions producing the power in the batteries generate heat [START_REF] Sato | Thermal behavior analysis of lithium-ion batteries for electric and hybrid vehicles[END_REF]. This heat generation increases as the electrical currents going through the batteries present higher values. The high temperatures can damage the batteries by causing fire hazards, or with less gravity, reducing their lifetime or performances [START_REF] Vetter | Ageing mechanisms in lithium-ion batteries[END_REF][START_REF] Allart | Gestion et Modélisation électrochimiques des batteries Lithium-ion[END_REF]. The thermal management of the batteries during the whole vehicle operation is then a key element to reach the demands of industrials and reach the performances allowing long and fast travels.

I.1.2 Battery thermal management systems

Tackling the thermal issues related to the Li-ion batteries operation requires a lot of effort to develop systems capable of cooling the batteries during the car operation. These systems are called Battery Thermal Management Systems (BTMS). Several technologies coexist in the EV market: phase change materials, indirect cooling, direct or immersion cooling (IC) or finally, two-phased liquid cooling. The role of the BTMS is threefold during the vehicle operation: keep the batteries in a safe range of temperatures, ensure a good temperature homogeneity between all the batteries in the pack and inside each battery. So far, the immersion cooling technology seems to be the most effective solution in terms of heat transfer performances [START_REF] Ehsan | Analyzing Thermal Management Methods of Li-ion Battery Modules[END_REF]. Also, the number of patents related to immersion cooling techniques rose significantly in the past years, as illustrated in

I.2 Immersion cooling of Li-ion batteries: state of the art and challenges

I.2.1 Challenges in battery thermal management systems

Since we have presented the main constraints that the electric vehicles and BTMS must fulfill, we review here the challenges induced by the physics occurring in such systems. Those systems feature complex physical behaviours during a normal operation of the vehicles.

First, the heat generated by the batteries during the operation of the vehicles, constituted by the alternation of charging and discharging sequences, can lead to the overheating of the batteries. The overheating hazard is more likely to occur as the battery pack is submitted to high electric currents. The fast charging and operation at high speeds of the cars typically require such values of the current. Then, with that usage, the batteries in the pack do not have much time to decrease their temperature.

As stated above, the batteries can reach high temperatures due to the heavy electrical currents submitted to the pack. By contrast, the cooling fluid entering the battery pack at the ambient temperature can provoke thermal shocks when the outside temperature is significantly lower. Also, the cooling fluid flow can present complex patterns due to packs arrangement geometry and the small cell spacing. In this perspective, the focus can be on the pressure of the cooling fluid inside the module by limiting the power consumption of pumps generating the flow. Then, the heat transfer can also be hard to assess because of the complex flow patterns and high-temperature gradients between the domains. Furthermore, more extreme events can occur in the BTMS. In some cases, the overheating can get so intense that it provokes thermal runaway [START_REF] Liu | Thermal issues about Li-ion batteries and recent progress in battery thermal management systems: A review[END_REF]. In such an event, one of the pack's battery reaches such a high temperature that it generates smoke and fire. As numerous cells are packed together, the heat spreads quickly through all the packs, igniting the other cells. This phenomenon features high-temperature gradients and occurs in a short time interval.

Finally, some BTMS use immersion cooling technology with a two-phased liquid-gas coolant flowing in the battery pack [START_REF] Wang | Thermal performance predictions for an HFE-7000 direct flow boiling cooled battery thermal management system for electric vehicles[END_REF]. This technique implies dealing with the pack's local phase changes and induced pressure variations.

A good representation of these physical processes is crucial to ensure a good design of the cooling systems. The numerical simulation stands as a powerful tool to tackle these challenges.

I.2.2 Simulation of immersion cooling systems for Li-ion batteries

Performing the simulation of immersion cooling systems, from the internal behavior of the Li-ion batteries to the heat transfer at a macroscopic scale, requires in practice to address a multi-physics problem. Several modeling approaches exist to perform the simulation of BTMS. The authors in [START_REF] Shahjalal | A review of thermal management for Li-ion batteries: Prospects, challenges, and issues[END_REF] identify the modeling methods to represent the thermal and electrical behavior of the Li-ion batteries in immersion cooling systems. They state that the usual approach is to consider the batteries and cooling system in either a 1D, 2D, 3D geometry or use a so-called lumped model. Then, different types of models can be found in the literature to represent the thermal and electrical behavior of the batteries.

For instance, the authors in [START_REF] Jaguemont | 1D-Thermal Analysis and Electro-Thermal Modeling of Prismatic-Shape LTO and NMC Batteries[END_REF] developed a 1D electro-thermal model to couple the electrical parameters of the battery, such as voltage and resistance, among others, with the temperature prediction of heated Li-ion batteries. The lumped body approach consider the temperature uniformly distributed in the battery. For example, [START_REF] Amini | A thermal model for Li-ion batteries operating under dynamic conditions[END_REF] apply a lumped thermal model to compute the heat generation in the battery and heat transfer between the battery surface and surrounding fluid. The lumped and 1D model developed in these studies featured a good representation of the temperature evolution with respect to experimental data.

On the other hand, the most common approach in the literature is to use Computational Fluid Dynamics (CFD) solvers to simulate the immersion cooling systems. At the scale of the battery pack, the heat transfer between the heated batteries and the cooling fluid, being either liquids or gas, is usually computed using Conjugate Heat Transfer (CHT). For most studies performing this approach, the Navier-Stokes equations are solved in the fluid domain. The heat equation is solved in the solid domain representing the batteries [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF][START_REF] Lin | A review on recent progress, challenges and perspective of battery thermal management system[END_REF]. These equations are coupled through temperature and heat flux continuity at the interface between both domains. The unknowns of interest are the temperature fields in fluid and solid domains, the fluid flow and pressure in the whole setup. The primary outcomes of studies performing the numerical simulation of immersion cooling systems can be summed up as follows: compute the temperature evolution at specific locations of the battery pack for given charging or discharging conditions, estimate the heat transfer coefficients to compare the cooling performances of some configurations depending on the type of coolant or the geometry of the arrangement of the cells.

The CFD studies in the literature are diverse and focus on different physical aspects of the immersion cooling systems. The numerical investigations of immersion cooling systems with CFD tools usually fulfill one of the two following objectives. The first objective consists of using the CFD models to characterize and analyze immersion cooling systems' typical heat transfer features in easy test case configurations. The second objective is to employ this numerical tool to directly assess the performances of some innovative cooling method or device from a pure design perspective.

For the first objective, the studies usually focus on the physical features listed in the following. Geometries of the arrangement of the cells in the modules are expected to affect the flow of the cooling fluid and heat transfer performances. The work performed in [START_REF] Wang | Thermal investigation of lithium-ion battery module with different cell arrangement structures and forced air-cooling strategies[END_REF] consists of a comprehensive investigation of this feature on the cooling efficiency. They assess the impact of the flow direction with respect to cylindrical cells set up in aligned, staggered or bunched geometries. In the same spirit, the study from [START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF] consists in finding the location of hot spots generated by their selected geometry in the system cooled by air. With small configurations, it is also easier to study the impact of the flow regime on the cooling. The authors in [START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF] evaluate the temperature evolution, experimentally and with a 2D CFD model, at specific locations of the cell, in natural and forced convection regimes. For instance, in work from [START_REF] Camilleri | Prediction of the Heat Transfer Coefficient in Direct Oil Cooling of Lithium-Ion Batteries[END_REF], they use CFD simulations of a bunch of cylindrical cells cooled with a high heat transfer coefficient liquid and compute the temperature for flow regimes featuring increasing Reynolds numbers. Finally, the type of fluid is also of importance. The difference is often made between air or liquid cooling systems in the literature. For further comparison, [START_REF] Chalise | Conjugate Heat Transfer Analysis of Thermal Management of a Li-Ion Battery Pack[END_REF] uses a small and straightforward geometry of the cooling system to perform CFD simulations with two fluids and assess the difference directly in terms of temperature resulting from the heating. The purpose of the above works is to represent heat transfer patterns in a small configuration and guess how the heat transfer would behave in larger configurations.

The second objective, on the opposite, considers more realistic immersion cooling devices represented numerically to assess the cooling efficiency of the whole set-up. Here, the CFD is used for the direct design of the systems. In this kind of studies, the geometry of the set-up is usually more complex or, at least, involves the representation of more numerous Li-ion cells.

The authors in [START_REF] Mahamud | Reciprocating air flow for Li-ion battery thermal management to improve temperature uniformity[END_REF] assess the performances of their reciprocating flow cooling system thanks to CFD simulation of the conjugate heat transfer. In this study, a complete battery pack is represented, even if some simplifications and symmetries alleviate the computational costs. The study's outcome stands in the assessment of the overall configuration thermal performances.

In [START_REF] Huat | Computational fluid dynamic and thermal analysis of Lithiumion battery pack with air cooling[END_REF], the authors built a whole battery module. The originality of their simulation approach is that they represent the full geometry, including the structural components of the battery pack. Their computation gives the temperature field for all those elements and provides a good overview of the system's behavior for engineers addressing its practical design.

An other example in [START_REF] Suresh Patil | A novel dielectric fluid immersion cooling technology for Li-ion battery thermal management[END_REF] shows an innovative method, mixing immersion and tab cooling (see section II.3.2 in chapter II for more details on cooling techniques). Even if the geometry comprises fewer cells, the CFD tool allowed the authors to test different configurations of their cooling device and assess the performances directly by evaluations of the temperature fields. Also, in [START_REF] Behi | A new concept of thermal management system in Li-ion battery using air cooling and heat pipe for electric vehicles[END_REF] the authors try a new type of cooling module mixing air and indirect tab cooling. The thermal performances are investigated experimentally. The 3D CFD model developed to represent the battery module is also used to perform further investigations by virtually changing the experimental geometry on cell spacing, ambient temperature or inlet air velocity on the heat transfer performances. This kind of studies presents a good merge of the two objectives stated above.

The review of these studies highlights some challenges usually encountered in the simulation of immersion cooling systems. First, as said earlier, real-life battery packs are composed of thousands of cells. Few studies address the simulation for such large geometries in immersion cooling configurations despite the objectives cited above. Representing the heat transfer in a real battery pack, at full scale with the thousands of cells and the complex induced structure, is still a challenge. The simulation of comprehensive configurations with precision, including the internal electrochemical process in the batteries, can be prohibitive with CFD models. Also, the meshing for these simulations can be challenging to realize. Industrials seek to set as many battery cells as possible in a given volume, making the spacing between them often tiny compared to the characteristic size of these devices, sometimes lower than the millimeter. The resolution of flow equations can be very challenging in these areas. Also, many of these studies used air as the cooling liquid. The coolant is usually a complex fluid with non-trivial thermal behavior in industrial applications. Also, the heat produced by the batteries comes from the internal chemical reactions occurring inside. So, representing the whole physics from the chemical level up to the module level's heat transfer would require setting up a multi-physics model coupling several levels. Only few studies, like for instance in [START_REF] Suresh Patil | A novel dielectric fluid immersion cooling technology for Li-ion battery thermal management[END_REF][START_REF] Jithin | Numerical analysis of single-phase liquid immersion cooling for lithium-ion battery thermal management using different dielectric fluids[END_REF], performed such a comprehensive approach.

Furthermore, the simulation of the immersion cooling systems with a two-phase fluid (liquid boiling) remains a difficult task. To all the constraints cited above, the model should represent the local phase changes in the whole geometry by capturing the moving boundary between the two phases.

Finally, the issues stated above can also explain the difficulty to find studies with a complete description of the multi-physics problem. The reproduction of results from literature can sometimes suffer from the lack of comprehensive information for the various parameters involved in such problems.

I.2.3 Uncertainty quantification for immersion cooling of Li-ion batteries

The immersion cooling systems come with a lot of uncertain parameters to consider. First, according to the scale considered to represent the heating of the batteries, some laws have to be set up for some key parameters governing the batteries' electrical and thermal behavior.

More specifically, the authors in [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF] identified the main scales that can be considered to represent these systems, from the batteries to the cooling system as a whole. First, the internal behavior of the batteries can be solved at the electrochemical level by solving electric charge conservation equations and the resulting heat production. The electrical behavior can also be computed at the macroscopic level of the cell by representing the electrical features of the cells. The heat generated is computed using phenomenological laws (such as Joule heating etc...) involving those electrical parameters [START_REF] Bernardi | A General Energy Balance for Battery Systems[END_REF]. Finally, the heat transfer between the batteries and the surrounding fluid can be solved for one cell or at the level of a whole module of batteries, involving many cells and complex flow resolution. Depending on which scale the models start from, assumptions have to be made to represent the physical behavior which would not be solved directly. All of these approaches consider models more or less distant from the internal physics of the cells. They are either constructed from experimental observations or physical considerations. Most importantly, they are defined through parameters coming with some epistemic uncertainty [START_REF] Der | Aleatory or epistemic? Does it matter?[END_REF], meaning that it is hard to choose values for these coefficients with strong confidence if no further information is provided.

These unknowns are expected to substantially impact the outcoming quantities of interest (QOI) when assessing the performances of the cooling systems. From the simulation perspective, these uncertainties should impact strongly the values predicted by the models. The predictive character of the numerical codes can be made more robust by dealing with those uncertainties. The forecasted quantities can lead to more insightful analysis and learnings when all these variabilities can be quantified.

An Uncertainty Quantification (UQ) approach can address these issues in numerical simulation problems. In this framework, simulation outcomes are not seen as a single deterministic output, but they give predictions in a given interval with an associated probability, which one aims to quantify. The accuracy of the models is assessed in terms of distributions and statistics (mean, variance, quantiles...) on these outputs. Specific mathematical tools need to be implemented and coupled to the numerical solvers to perform such an approach.

Once one can represent the uncertainties of the input parameters, some standard methods are performed in literature to quantify the uncertainties on the numerical model's responses.

The statistics resulting from the uncertainties can be computed using numerous samples of the inputs. This is called the Monte Carlo (MC) method adapted to computer models [START_REF] Caflisch | Monte Carlo and quasi-Monte Carlo methods[END_REF][START_REF] Kroese | Handbook for Monte Carlo methods[END_REF]. Each sample corresponding to an evaluation of the solver, surrogate models, are usually employed to alleviate the computational costs associated with numerical models, which are often expensive in practice [START_REF] Sudret | Surrogate models for uncertainty quantification: An overview[END_REF]. Further analysis on the model's behavior can be assessed with methods such as Global Sensitivity Analysis (GSA) [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF]. Here, the impact of the uncertain inputs on the variability of quantities of interest can be estimated [START_REF] Razaaly | Impact of geometric, operational, and model uncertainties on the non-ideal flow through a supersonic ORC turbine cascade[END_REF]. Finally, an inverse problem can be solved to infer some uncertain inputs, using the information provided by some trustful data, such as experimental measurements. This way, the models can be calibrated under a probabilistic framework, following the Bayes paradigm, and the prior uncertainties on the inputs are reduced. This method can be used to reconstruct values and distributions of inputs, using the knowledge brought by the observation of some output data [START_REF] Kennedy | Bayesian calibration of computer models[END_REF][START_REF] Francesco | Predictive numerical simulations for rebuilding freestream conditions in atmospheric entry flows[END_REF]. Some references have employed uncertainty quantification approaches in the simulation of Li-ion batteries. Most of these authors focus on the equations governing the electrochemical behavior within a Li-ion cell [START_REF] Vahid Esfahanian | Capacity Fade Analysis of Lithium-Ion Cells Utilizing Uncertainty Quantification Approach[END_REF][START_REF] Hadigol | On Uncertainty Quantification of Lithium-ion Batteries: Application to an LiC 6 LiCoO 2 cell[END_REF][START_REF] Laue | Model-Based Uncertainty Quantification for the Product properties of Lithium-Ion Batteries[END_REF][START_REF] Andrew | Stochastic capacity loss and remaining useful life models for lithium-ion batteries in plug-in hybrid electric vehicles[END_REF]. Some authors [START_REF] Tong | Correlating uncertainties of a lithium-ion battery -A Monte Carlo simulation: A Monte Carlo simulation of a lithium-ion battery model[END_REF][START_REF] Chen | Modeling and characterization of the mass transfer and thermal mechanics of the power lithium manganate battery under charging process[END_REF] focus on coupling electronics partial differential equations to the thermal equation to compute the heat produced within a single Li-ion cell. Not many references are studying the effects of uncertainties for the heat transfer between the batteries and cooling fluid in immersion cooling systems.

I.3 Motivations and contributions of the thesis

I.3.1 Industrial collaboration with Exoes

In this industrial context, this thesis project has emerged from the collaboration with the company Exoes in Gradignan, France [START_REF] Exoes | Avenue de la Grande Lande[END_REF]. Exoes develops battery modules with immersion cooling systems for car manufacturers. Exoes provides a battery pack solution responding to customer needs in terms of space available in the car, power required by the engine, and cost limitations. They propose the suitable shape and type of Li-ion batteries and the integrated immersion cooling system. Then, their expertise is about the construction and the analysis of these cooling systems.

The collaboration with Inria comes naturally from the need of Exoes to develop advanced numerical tools to perform more profound analysis of these systems and support the design of new systems. The complexity of the immersion cooling systems requires the development of high-fidelity solvers dedicated to the simulation of these systems and rigorous treatment of uncertainties, which are the main focus of this thesis. Moreover, the experimental data generation from Exoes allowed a systematic validation of the numerical models.

I.3.2 Thesis motivations

Considering the aspects described above, from the industrial design of immersion cooling systems to their simulation and the treatment of the uncertainties, the objectives of this thesis project can be described as follows.

Numerical simulation of immersion cooling systems for Li-ion batteries

The first objective of the thesis is to develop numerical models to perform the simulation of immersion cooling systems for Li-ion batteries. Models should predict thermal quantities, such as the transient evolution of temperature under different conditions. A trade-off between computational cost and fidelity of the models will be explored: two models of different levels of fidelity will be considered. They will be designated in the following, respectively, as low fidelity (LF) and high fidelity (HF) solvers. The nature of the results produced by these models can be put in perspective with the respective numerical resources involved.

Identify and characterize the uncertainties coming from immersion cooling systems for

Li-ion batteries The multi-physics system requires dealing with uncertainties to perform informative and predictive simulations. This objective consists of producing a physical analysis of the immersion cooling systems to identify which parameters should be considered uncertain. Based on literature reviews and experts analysis from Exoes, the list of such parameters can be established. The collaboration with Exoes allows bringing much information for this particular task and is essential to characterize precisely system's uncertainties. The uncertainties to be considered (including range of variation and distribution) should be defined for the selected parameters at this point.

Robust validation and calibration of the models with experimental data

Once the numerical models are set up, and the uncertainties are well characterized, the final objective is to produce simulation results to get analysis elements and design the cooling system, including more profound knowledge of its behavior. First, the numerical models should be validated against experimental data. Using experimental data directly can permit calibrating the model under a probabilistic framework (Bayesian calibration). The outcomes of such analysis will lead to informative values of the uncertain parameters, allowing model responses close to the experimental data. In this way, uncertainties in the numerical prediction can also be reduced. Overall, the design of the immersion cooling can be addressed with robust analysis elements, merging information provided from trustful experimental data and numerical predictions quantifying the level of the unknown in the model's responses.

I.3.3 Scientific contributions of the thesis

According to the objectives mentioned just above, the contributions of the PhD project to tackle these issues can be listed under the four following points.

Fast and low fidelity solver adapted from Exoes' solution First, a fast numerical model, called ICExo, developed originally at Exoes, was modified to simulate batteries immersed in a cooling fluid. The numerical model allows the prediction of temperature and electrical quantities such as voltage and state of charge of the batteries under varying charging and discharging conditions. This solver is seen as a low fidelity model as it does not solve the complete set of thermal and flow equations as it would be done in a pure CFD approach.

The contribution of this solver stands in its low computational cost and still good predictive character. From a practical point of view, the first version of the model developed at Exoes was not adapted for a suitable usage in an uncertainty quantification framework and to be run with multiple input conditions in parallel on the cluster PlaFRIM [START_REF]PlaFRIM experimental testbed, supported by Inria[END_REF]. A part of the work was to start from the initial version of the code and adapt it to reduce the computational costs. In the end, the code was able to run with fewer resources, featuring a reduced amount of required data and an optimization of the elementary operations.

High fidelity open-source CFD solver applied to immersion cooling systems Secondly,

an open-source CFD model, TrioCFD [140], developed at CEA Saclay, France, was selected to be used for the first time in battery cooling applications. Based on the Finite Element Method (FEM), the model showed good capabilities to simulate immersion cooling systems.

The numerical solver was used to solve the coupled problem, constituted by the transient conjugate heat transfer between the heated solid batteries and the cooling fluid. As this code was used for the first time for this application, a significant part of the work consisted in selecting immersion cooling test cases and adapting the CFD code to the references' setups.

Usually, the nuclear applications present significant differences in the order of magnitudes for the physical processes compared to Li-ion batteries immersion cooling problems. Then the reproduction of results from the literature represented a large part of the investigations performed during the thesis. The model's accuracy was demonstrated by comparing it with commercial software and experimental data. Uncertainty quantification methods such as surrogate-based construction and global sensitivity analysis were performed. Finally, Bayesian calibration of the uncertain input parameters was addressed to validate the numerical model with respect to experimental measurements conducted directly by Exoes. Thanks to this process, the uncertainty in the inputs values and the temperature prediction are reduced significantly. Overall, the numerical model and the uncertainty approach gave some informative insights on the immersion cooling behavior and the values of the involved physical parameters.

Uncertainty propagation and calibration with low fidelity solver

Internal resistance model accuracy assessment based on CFD simulations and experimental measurements

The last contribution of the work is about an uncertainty quantification approach performed using the data provided by the high fidelity CFD solver. The focus is set here on the internal resistance of the batteries, being one of the most important parameters when modeling the heat produced by Li-ion cells. The test case addressed consists of the heating of batteries, submitted to a constant electrical current and cooled with the surrounding air. While solving the conjugate heat transfer problem, a model is constructed to represent the behavior of the internal resistance, respecting physical constraints taken from literature. The constructed model allowed to deal practically with the uncertainties inherent to the resistance parameter. The forward propagation of the uncertainties allows assessing the impact of the input uncertainties on the CFD model's temperature prediction. Finally, a Bayesian calibration is performed using experimental data from the test case. Values of the internal parameters leading to a numerical response close to the experimental data are learned from this process. The uncertainty in the temperature prediction was also considerably reduced in that configuration. From a practical point of view, this overall approach constituted a method to reconstruct the behavior of the resistance from measurements of the temperature of the heated Li-ion batteries.

I.4 Outline of the manuscript

The present manuscript is structured as follows.

Chapter II introduces the Li-ion batteries operation, the parameters used to characterize their behavior and the main physical processes at stake. Thermal issues related to batteries' operation are also covered, highlighting the need of using BTMS in Electric vehicles applications. The different types of BTMS are described, emphasizing the system of interest in this work, the immersion cooling technique. The modeling approaches existing in literature to simulate the batteries and their heat production are also covered. The usual approach for the simulation of conjugate heat transfer in this configuration is also briefly described at this point.

Chapter III is devoted to the description of the LF solver from Exoes, called ICExo. The full operation of the model is provided, including the simulation of thermal and electrical features. An original experimental test case conducted at Exoes is presented. The data outcoming of this experimental case is used to validate the behavior of the numerical model under different conditions.

Chapter IV consists of the presentation of the CFD solver. The equations solved in the framework of this thesis are fully described and given according to the corresponding physical justifications. The code behavior is first validated with a comparison against FLUENT on a test case of forced convection. Then the code is further validated against experimental data. This study highlights the importance of the source term selection and the impact of 2D or 3D computations in the immersed batteries temperature prediction submitted to a constant discharging current.

Chapter V covers the theoretical framework to address uncertainty quantification problems. First, the overall point of view of uncertainty quantification approaches is given. The mathematical tools required in the usual UQ framework are briefly described. Then we focus on the description of the methods used in the applied UQ studies performed in this work.

Chapter VI details the steps and the results of the contribution mentioned above about the uncertainty quantification study using the data provided by the LF model. The main results of this chapter are the calibration of the uncertain multi-physics parameters, the resulting uncertainty propagation and the outcomes of the Global Sensitivity Analysis performed in a realistic racing case.

Chapter VII describes the last contribution about the internal resistance model accuracy assessment. After briefly recalling the test case, the process to construct the internal model is given. Next, the Bayesian calibration shows the uncertainty reduction capabilities of this approach from the propagation of prior and then posterior distributions of the resistance model parameters. The overall process shown in the chapter constitutes the method for reconstructing the internal resistance model from the measurements of temperature on the Li-ion cells surface. Finally, an exploratory numerical test is performed considering the effect of the state of charge (SOC) in the internal resistance model, added up with the temperature dependency. The proposed method allows assessing the impact of the low state of charge on the temperature evolution.

Chapter VIII provides some concluding remarks on the results of the thesis. Perspectives for short and long term actions are given to draw a roadmap of the work required to go deeper in analyzing the current results and explore new directions.

LITHIUM-ION BATTERIES AND IMMERSION COOLING

Chapter abstract

This chapter covers the physics of Li-ion batteries and cooling technologies. First, the different components and characteristics of batteries are reviewed. The electrochemical phenomena governing the operation of Li-ion batteries are described from a physical perspective. Particular attention is paid to the thermal issues related to the operation of the Li-ion cells. A review of the models used in literature to represent these systems' electrical and thermal behavior is presented. Secondly, the existing cooling techniques for battery thermal management are introduced. Physical explanations of the processes attempt to illustrate the good performances offered by immersion cooling. The physics of the immersion cooling technology is detailed, and the equations used to solve the conjugate heat transfer problem are expressed. High specific power means that the car can make strong accelerations and reach high speeds.

Outline

These batteries are capable of delivering and receiving high electric currents. The high specific energy is traduced by a great range of autonomy for the vehicles. The following diagram highlights the good performances of Li-ion batteries. The red beam concerning Li-ion batteries lives between high specific energy and high specific power range. This technology is the closest to the capabilities currently offered by internal combustion engines (IC Engine in the figure).

They also present other advantages than the energy and power density, which makes them practical to use: a great shelf life, with a relatively low rate of capacity loss over time.

Also, the temperature range of these batteries is extensive, usually from -40 o C to 70 o C , making them a practical choice when designing vehicles for diversified applications.

The following paragraphs detail the main components of these systems. A Li-ion battery is composed of: a positive and negative electrodes, a separator, an electrolyte [START_REF] Allart | Gestion et Modélisation électrochimiques des batteries Lithium-ion[END_REF][START_REF] Zhengming | Lithium-Ion Battery Systems and Technology[END_REF]. These components are disposed in a layered structure under various existing geometries, which are detailed later within this manuscript.

Negative electrode

The active material of this electrode is often made of graphite carbon.

The material composing the anode should present physical features as a minimal volume expansion and stress. Also, it should be chemically stable under an extensive range of temperatures during the charging and discharging sequences of the battery [START_REF] Zhengming | Lithium-Ion Battery Systems and Technology[END_REF]. The chemical structure of the graphite material allows the insertion of the Lithium during the charging process [START_REF] Allart | Gestion et Modélisation électrochimiques des batteries Lithium-ion[END_REF]. It is associated with a copper (Cu) plate to collect the electric current. The active material hosting the Li elements is dropped on the copper collector.

The surface of the electrodes are covered with a film called SEI for Solid Electrolyte Interface [START_REF] Vetter | Ageing mechanisms in lithium-ion batteries[END_REF][START_REF] Kurzweil | Lithium Battery Energy Storage[END_REF]. This last one is the anode during the discharge sequence.

Positive electrode

The active material is usually metal oxide. The electrical current collector is made of aluminium (Al). The active material capable of letting the ions Li + to insert in is dropped on this collector. This one is the cathode during the discharge sequence.

Note that the choice of materials for cathode and anode will have a direct impact on the energy storage capacity of the cell, depending on the compactness of their chemical structure [START_REF] Allart | Gestion et Modélisation électrochimiques des batteries Lithium-ion[END_REF].

Electrolyte

The electrolyte is the aqueous solution allowing the ions Li + to flow from an electrode to each other. It is carefully chosen depending on the materials of the electrodes as it should not provoke chemical reactions aside from the reduction-oxidation. The most commonly used salt materials are LiPF 6 , LiBF 4 . The electrolyte should also present a high dielectric constant [START_REF] Zhengming | Lithium-Ion Battery Systems and Technology[END_REF].

Separator The separator is a porous membrane disposed between the electrodes. Its role is to prevent the electric contact between the electrodes and allow the ionic flow in the meantime. The microporous membranes are usually made of polyethylene or polypropylene, or both.

II.1.2 Electrochemical reactions

As we just detailed the role and chemistry of the structural elements composing a Li-ion battery, let's review the chemical equations on both electrodes. The operation of a Li-ion battery is governed by the oxidation-reduction equations, generating the electric current and flow of Li + ions in the electrolyte. For the discharge, the half-chemical reactions occurring between the graphite carbon C from negative electrode, and the metal oxide Li x M O 2 of the positive electrode reads [START_REF] Allart | Gestion et Modélisation électrochimiques des batteries Lithium-ion[END_REF]:

DISCHARGE REACTION (II.1) anode -deintercalation Li x C 6 → 6C + xLi + + xe - cathode -intercalation Li 1-x M O 2 + xLi + + xe -→ LiM O 2 (II.2)
A schematic representation of these electro-chemical equations is given in

II.1.3 Electrical characterization

It is possible and sometimes more relevant to see the Li-ion batteries from a broader perspective than the electro-chemical scale. Then the batteries can be fully characterized by parameters describing the electrical state of the whole chemical structure.

Open circuit voltage (OCV)

This one is the nominal potential of the cell, expressed in volts, often denoted as OCV. It represents the electrical potential difference between the positive and negative electrodes when nothing is connected to the battery. It is the equilibrium voltage of the cell at rest. Each Li-ion cell is characterized by its open-circuit voltage. This parameter is directly linked to the state of charge of the battery, following non-trivial relationship to assess in practice [START_REF] Zhang | A Generalized SOC-OCV Model for Lithium-Ion Batteries and the SOC Estimation for LNMCO Battery[END_REF][START_REF] Yu | A Comparative Study on Open Circuit Voltage Models for Lithiumion Batteries[END_REF].

State of charge (SOC)

Usually denoted SOC, the state of charge indicates how much the battery is still capable of delivering electrical current. It is a unitless number, often expressed in percentage. A SOC equal to 100% means the battery is fully charged. When the SOC is at 0, the battery is said empty. As said above, the OCV and SOC are related: during the battery operation, when the OCV is at its maximal value, the SOC is 100%, and the battery is full.

On the opposite, when the voltage of the battery drops to its minimal value, the SOC is 0%, and the battery is fully discharged. Many research works investigate the link between these two quantities and then address the issue of giving a reliable method to compute the state of charge [START_REF] Chang | The State of Charge Estimating Methods for Battery: A Review[END_REF]. However, an agreement on the typical relationship between the OCV and the SOC seems to be established in the literature, similar to the one provided by [START_REF] Tingting | Analysis on the influence of measurement error on state of charge estimation of LiFePO4 power Battery[END_REF] and

illustrated in Fig. II.4.

State of health (SOH)

This parameter is used to quantify the age of the battery. First, the batteries present deprecated performances over time, even if not used. Also, Li-ion batteries are given to operate for a fixed number of charge/discharge cycles throughout their life.

The capacity of the batteries decreases with time and the number of performed cycles (see

Fig. II.5
). Different chemicals phenomena occurring within the batteries are causing this loss [START_REF] Kurzweil | Lithium Battery Energy Storage[END_REF]. Some molecules composing the solvent intercalate into the electrode materials, preventing the Li + ions from being stored in the molecular structure. The surface film on the electrodes may be modified, consuming some Lithium ions. The state of health is then used to consider how much battery capacity has been lost over time compared to day one.

Discharge and charge rate (C-rate)

This unitless quantity, commonly written C-rate, is used when regarding the operation of the battery and does not characterize its state. Depending on its original capacity C, it gives the intensity of the electrical current used to charge or discharge a battery. This quantity is helpful to compare under which intensity batteries are used, even if they are not of the same type and present different characteristics. For instance, discharging a battery at a 2 C-rate means that the submitted electric current is such that the full discharge of the battery will occur in 30 minutes. Charging a battery at 0.5 C-rate means that a full discharged battery would be filled up in 2 hours. The capacity C of the battery is then stated as the 1 C-rate. The link between the time of charge/discharge expressed in hours, here denoted τ , and the C-rate C r is:

τ = 1 C r (II.3)
Knowing the capacity C of the battery, the applied C-rate C r allows to retrieve the associated electrical current I using:

I = C r • C (II.4)
This parameter is often used to compare the evolution of some quantities of interest related to the battery characterization (voltage, capacity,...) under different discharging conditions, as illustrated in Fig. II.6.

Internal resistance (R)

The layered physical elements composing the cell structure are crossed by electrons during the cell operation. Each of these presents an inherent electric resistivity and obstruct the electrons transfer. As for every ohmic electrical device, the resistance directly relates to the voltage. Due to the composition and the chemistry of the cell, the internal resistance depends on parameters such as the internal temperature, the state of health, the state of charge and the employed C-rate during operation. The internal resistance also depends on the electric current submitted to the batteries. These parameters govern many phenomena occurring during a battery operation, as it is responsible for voltage drop and energy loss.

II.1.4 Types of Li-ion batteries

There is a wide range of applications where Li-ion batteries are used; many geometries and types of batteries coexist in the market. The differentiation is based on the geometries and the materials used for the positive electrode. The most common geometries used in the automotive industry are: prismatic, pouch, and cylindrical cells [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF]. A specific geometry is used depending on the design of the electric vehicle under development. Each of these geometries, visible in Fig. II.7, presents advantages and drawbacks, primarily related to thermal management issues.

The difference between batteries are also characterized by the material used for the positive electrodes, as most of them are made with a graphite negative electrode. In [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF] As highlighted in Tab. II.1, most of the battery types present similar voltage and capacities. Other more practical characteristics are then considered when choosing the materials, such as the power delivered by the cell, safety, cost and lifetime. Lithium Iron phosphate LiF eP O 4 seems to show the best behavior considering these aspects [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF]. This chemistry also presents characteristics slightly different than the others as a significant hysteresis and a lower open-circuit voltage range, as stated in [START_REF] Ceraolo | Experimentally-Determined Models for High-Power Lithium Batteries[END_REF]. The advantages and drawbacks of each chemistry are reviewed in detail in [START_REF] Kurzweil | Lithium Battery Energy Storage[END_REF].

II.1.5 Thermal issues

II.1.5.1 Heat generation processes

The heat is generated by the electrochemical reactions between all the different layered elements composing the Li-ion cells. Several phenomena are mainly responsible for the heat generation in a Li-ion battery, considering the usual time scales of the study. There are two types of heat generation occurring in the battery: reversible heat and irreversible heat.

First, the reversible heat is caused by the entropy changes coming from the electrochemi-cal reactions. This heat accounts for the absorbed energy when the lithium ions intercalate in the negative electrode during charge. This reaction is endothermic during charge. During the discharge, the opposite reaction occurs and is exothermic [START_REF] Kang | Thermal Analysis of a Parallel-Configured Battery Pack (1S18P) Using 21700 Cells for a Battery-Powered Train[END_REF]. The exact amount of energy is released or absorbed by the battery, during discharge or charge respectively [START_REF] Sato | Thermal behavior analysis of lithium-ion batteries for electric and hybrid vehicles[END_REF].

The irreversible heat sources are coming from the ohmic heating and polarization reactions. First, the ohmic heat source is generated by the internal ohmic resistance of the battery.

It's caused by the ions and electrons transport across the solid and electrolyte phases [START_REF] Du | An investigation of irreversible heat generation in lithium ion batteries based on a thermo-electrochemical coupling method[END_REF].

As it is explained in [START_REF] Sato | Thermal behavior analysis of lithium-ion batteries for electric and hybrid vehicles[END_REF], this heat generation is dependent on the internal resistance of the battery. And the overall internal resistance is caused by the materials composing the layered structure of the battery: the current collectors, the active electrode materials, the separator and the electrolyte. When submitted to an input current, the transfer of charged particles in the battery is obstructed by these materials and the energy losses are released under heat form. This phenomenon is the Joule heating.

Secondly, the polarization heat generation is caused by energy losses during the irreversible electrochemical reactions [START_REF] Du | An investigation of irreversible heat generation in lithium ion batteries based on a thermo-electrochemical coupling method[END_REF]. Those reactions are irreversible as they represent the degradation of electrodes surface due to the motion of ions [START_REF] Sato | Thermal behavior analysis of lithium-ion batteries for electric and hybrid vehicles[END_REF]. Then, this heat is directly related to the potential of the electrodes, as seen later in the electrochemical model of heat generation.

II.1.5.2 Effects of temperature

The temperature has a substantial impact on the batteries' performances. Indeed, all of the electrical parameters described above in the manuscript significantly depend on temperature. The temperature impacts two main features of the Li-ion batteries: ageing and thermal runaway.

The ageing is impacted at both high and low temperatures [START_REF] Vetter | Ageing mechanisms in lithium-ion batteries[END_REF][START_REF] Allart | Gestion et Modélisation électrochimiques des batteries Lithium-ion[END_REF]. Dendrites can form between the negative electrode and the electrolyte at very low temperatures. This fact might increase the resistance at the interface between those two media. The electrolyte presents a drop of ionic conductivity for such temperatures, slowing the electric charge transfers between the layered structure. On the opposite, at high temperatures, the solid electrolyte phase might melt within the electrolyte, and the negative electrode reacts again with the electrolyte. An additional heat source is created from these reactions. Also, the separator is susceptible to melt at high temperatures. The high temperatures also damage the structure of the electrodes.

The thermal runaway is reached when the cell's temperature is so high that heat dissipation is not possible anymore. It can have several causes, such as mechanical, electrical or thermal abuses [START_REF] Feng | Thermal runaway mechanism of lithium ion battery for electric vehicles: A review[END_REF]. This temperature is generally the fusion temperature of the Lithium elements. When such an event occurs, the gases released by the chemical reactions described above increase the pressure within the battery structure. Smokes are released through the safety exit valve together with the ejected electrolyte. Then, the battery explodes because of the highly reactive elements composing its structure. The triggered fire burns all the left materials.

The author in [START_REF] Allart | Gestion et Modélisation électrochimiques des batteries Lithium-ion[END_REF] proposes an informative summary of the temperature effects on the battery structure, summed here in Tab. II. This brief review of the various and undesirable thermal effects on the batteries' composition underlies the need for a thermal management system to control the batteries' temperature during operation. The battery thermal management systems are reviewed in II.3.

II.1.5.3 Parameters governing the temperature field evolution

Since we have reviewed the main physical mechanisms of heat generation within the Liion batteries, we can focus on the parameters governing the temperature distribution and evolution. They depend on the inherent chemistry and structure of the batteries. The two parameters governing the temperature evolution and distribution within the cells are the heat capacity, and the thermal conductivity [START_REF] Drake | Thermal conduction and heat generation phenomena in li-ion cells[END_REF].

The heat capacity will govern at which rate a volume is heated or cooled. The more significant the heat capacity, the more energy needs to be supplied to change the body's temperature. This parameter impacts the temperature evolution in transient heat transfer problems. The thermal conductivity rules the temperature distribution within a volume. It governs the temperature homogeneity of a medium subject to heat transfer. If the thermal conductivity is high, the temperature field tends to be more homogeneous, as the heat is well-spread thanks to the excellent conductivity of the material.

Due to the reactions occurring in the porous media of the battery, the thermal properties of the components are changing during operation and life of Li-ion cells. The proper way to evaluate those properties is to consider the 'effective' thermal conductivities and heat capacity, depending on the volume fraction of each chemical component living in the porous layers.

For the heat capacity, the authors from [START_REF] Drake | Measurement of anisotropic thermophysical properties of cylindrical Li-ion cells[END_REF] proposed an original approach mixing experimental and analytical methods to assess the value of this parameter in some Li-ion batteries.

They measured the heat capacity of a cylindrical battery, heated with a known source term, axially and radially with respect to the cylindrical geometry. The heat capacity in the axial test appeared to be lower than in the radial test. For cylindrical batteries of type 26650 (radius of R = 26 [mm] and height L = 65 [mm]) the radial test gave a heat capacity of

1605 J • kg -1 • K -1 and the axial test 1895 J • kg -1 • K -1 .
Finally, the heat capacity resulting from the axial test appeared to be the most accurate regarding the structural composition of the battery.

Furthermore, the thermal conductivity of Li-ion cells presents an anisotropic behaviour.

This property comes from the structural composition of electric Li-ion cells [START_REF] Saw | Electrochemical thermal analysis of 18650 Lithium Iron Phosphate cell[END_REF]. As seen earlier, all existing geometries of Li-ion cells (cylindrical, prismatic and pouch cells) present a layering of positive/negative electrodes, separator and electrolyte. As it can be understood from the figure II.2, the thermal conductivity of a Li-ion cell depends on the materials which compose the layering and their arrangement. It has been observed that the thermal conductivity is low in the direction crossing the layers and higher in the direction along a layer.

Many researchers has presented innovative methods, including experimental and numerical simulations, to evaluate this parameter [START_REF] Drake | Thermal conduction and heat generation phenomena in li-ion cells[END_REF][START_REF] Wei | A comprehensive study on thermal conductivity of the lithium ion battery[END_REF][START_REF] Huang | Experimental measurement of anisotropic thermal conductivity of 18650 lithium battery[END_REF][START_REF] Murashko | Thermal parameters determination of battery cells by local heat flux measurements[END_REF][START_REF] Loges | Thermal characterization of Li-ion cell electrodes by photothermal deflection spectroscopy[END_REF]. Regardless of the battery geometry considered, the conductivity in the cross direction of the layers is much lower than the one in the parallel direction of the layers. Indeed, these studies show that the cross layers thermal conductivity is around ten times smaller than the axial thermal conductivity. [START_REF] Wei | A comprehensive study on thermal conductivity of the lithium ion battery[END_REF] recently proposed a method to compute the thermal conductivities, considering the geometry of the layering and the volumetric fraction of the chemical components and highlighting differences between cross plane and axial conductivities. The values are sorted according to the choice of the negative electrode material and presented in Tab. II.3.

Battery type

Axial thermal cond. Cross layers thermal cond.

W • m • K -1 W • m • K -1

II.2 Modeling approaches for Li-ion batteries and heat source II.2.1 Modeling the characteristics of Li-ion batteries

There are several approaches to simulate the batteries' internal electric behavior and characteristics, depending on the scale considered. First, we detail the principles of modeling at an electro-chemical scale. Then, the approach with a larger scale is detailed to represent the main parameters involved in the electrical phenomena.

II.2.1.1 Electro-chemical scale

The behaviour of the battery can be represented at the chemical level by solving the partial differential equations governing the electrons and ions transfer within the layered structure of the battery [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF][START_REF] Basu | Coupled electrochemical thermal modelling of a novel Li-ion battery pack thermal management system[END_REF]. In [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF] a summary of the essential equations and variables at this scale is presented. This approach consists of a set of partial differential equations to solve, which represent the following physical processes:

• Electrochemical kinetics: the Butler-Volmer equation describes the interface between the electrodes particles (positive and negative) and the electrolyte [START_REF] Du | An investigation of irreversible heat generation in lithium ion batteries based on a thermo-electrochemical coupling method[END_REF]. It gives the

local current density ([A • m -2 ]
) function of the transfer coefficient for anodic and cathodic currents and the potentials of each electrode.

• Phase transition and ion transport -Solid-phase conservation of ions Li + : diffusion equation solving the concentration of Li + species in the solid phase.

-Electrolyte phase conservation of ions Li + : diffusion equation solving the concentration of Li + species in the electrolyte phase. A term is added to represent the transfer of Li + ions due to the flow of solvent and the transfer current, denoted j Li , coming from the insertion (or de-insertion) of the Li + ions in the electrode/electrolyte interface.

• Energy dissipation -Charge conservation in solid phase: Fick's law linking the transfer current and the electrical potential through the effective conductivity.

-Charge conservation in electrolyte phase: Fick's law of the ionic transport in electrolyte.

II.2.1.2 Equivalent circuit models

Another way, coming with a lower fidelity in the representation of the internal physics of the cell, is to represent the battery using an equivalent electrical circuit model. The battery operation is approximated with a suitable connection of classical electric components. This approach is preferred by electrical engineers who need simpler models and are more used to thinking in terms of voltage and current than chemical equations [START_REF] Ceraolo | Experimentally-Determined Models for High-Power Lithium Batteries[END_REF]. It gives the relationship between input and output parameters of interest related to the battery operation.

Modeling the behavior of the battery is then reduced to the representation of the electrical behavior of those components and their connections. Several electrical circuits are proposed to simulate the battery behavior based on empirical observations. Here, we review the most common networks used to represent the electrical behaviours of Li-ion batteries. The main features of these models are reviewed in [58] as follows:

II.2.2 Heat source term

The internal thermal behavior of the battery is represented through the heat equation solved in the domain, seen as a homogeneous or multi-phase medium. Depending on the modeling scale considered, the heat source term will result from different physical processes.

II.2.2.1 Heat generation at electro-chemical scale

At the electro-chemical scale (see II.2.1.1), the variables involved in the system of equations consider the transfer of species between the layers of the battery. The heat source term then depends on those variables and allows to compute the heat generation locally in the domain.

According to this approach, the heat source term is usually the result of three contributions [START_REF] Basu | Coupled electrochemical thermal modelling of a novel Li-ion battery pack thermal management system[END_REF][START_REF] Du | An investigation of irreversible heat generation in lithium ion batteries based on a thermo-electrochemical coupling method[END_REF]:

q = Q rev + Q p + Q ohm (II.5)
The expressions of these source terms involve the following unknowns and quantities: i s , i e are the current density in solid and electrolyte phases, respectively, i n is the transfer current connecting the electrolyte and solid phases, a s is the specific surface of the solid phase. Also, φ s , φ e and E eq are the potentials in the solid, electrolyte phases and at equilibrium. x is the vector containing the plan coordinates of the location considered within the domain. T is the temperature field in the battery, function of space coordinates and time.

As seen in II.1.5.1, Q p is the polarization heat generation due to over potential or the voltage loss caused by the degradation of electrodes surfaces [START_REF] Du | An investigation of irreversible heat generation in lithium ion batteries based on a thermo-electrochemical coupling method[END_REF].

Q p = a s i n (φ s -φ e -E eq ) (II.6)
Q rev stands for the reversible heat generation, due to the exothermic or endothermic reactions (charge or discharge) caused by the entropy variation.

Q rev = a s i n T ∂E eq ∂T (II.7)
Q ohm is the Ohmic heat generation due to the resistance to the electronic transport in the solid phase, and the resistance to ionic transport in the electrolyte phase.

Q ohm = -i e ∂φ el ∂x -i s ∂φ s ∂x (II.8)
This approach allows us to understand the heat generation process with its spatial dependency within the battery as the variables involved in the heat source are computed at each location of the domain. As exposed in Fig. II.9, one can observe the link between the current density, the heat generation and then the temperature distribution within the cell structure.

This method is seen as the highest-fidelity level, as it represents most of the electrochemical processes occurring in the battery and computes the outputs of interest with space and time dependency. The induced computational cost is then very high, and this approach might be prohibitive when one seeks to simulate a configuration with many batteries and surrounding fluid, as detailed later in this chapter.

Figure II.9: 3D fields of current densities, resulting volumetric heat and temperature distributions in the case of a cylindrical Li-ion cell. Figure courtesy of [START_REF] Basu | Coupled electrochemical thermal modelling of a novel Li-ion battery pack thermal management system[END_REF].

II.2.2.2 Heat generation at cell scale

Following the philosophy of the equivalent circuit model scale, the heat source can be computed directly using the electrical parameters' macroscopic values. The heat generation is seen at the level of the entire cell. With this point of view, the parameters involved are not computed from the electro-chemical reactions but seen as physical parameters playing a role in the whole solid cell. The fidelity of this approach is decreased as the microscopic phenomena are represented by more general terms in the heat generation equation. However, it turns out to be more practical when addressing the simulation of heated batteries from a larger perspective.

A model taking into account two main electro chemical phenomena was first proposed by [START_REF] Bernardi | A General Energy Balance for Battery Systems[END_REF] and widely used in thermal analysis of Li-ion batteries [START_REF] Cao | Thermal modeling of full-size-scale cylindrical battery pack cooled by channeled liquid flow[END_REF][START_REF] Xia | Thermal Analysis and Improvements of the Power Battery Pack with Liquid Cooling for Electric Vehicles[END_REF][START_REF] Karimi | Thermal management of lithium-ion batteries for electric vehicles[END_REF][START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF][START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF][START_REF] Damay | Thermal modeling and experimental validation of a large prismatic Li-ion battery[END_REF][START_REF] Kleiner | Thermal Modelling of a Prismatic Lithium-Ion Cell in a Battery Electric Vehicle Environment: Influences of the Experimental Validation Setup[END_REF].

According to this approach, the heat source is seen as the sum of two contributions. The first one is the heat generated by the entropy changes caused by the electrochemical reactions, denoted here as q e . This is the reversible heat described earlier in the manuscript.

The second contribution is the heat generated by the Joule effect, q J . It is caused by the electrons transfer across the layers composing the battery. Each layer presents an inherent electric resistivity and obstructs the electrons motion. The induced energy loss is then released in the form of heat. The expression from [START_REF] Bernardi | A General Energy Balance for Battery Systems[END_REF] in simplified form and commonly used reads:

q = q J + q e = I(E oc -U ) ohmic heating -IT dE oc dT entropic heating (II.9)
where E oc is the open circuit voltage of the battery, U its voltage, I the current load, T the temperature within the battery. Then, these terms can be computed using one of the equivalent circuit model. We will see in next chapters that the ohmic heating is usually computed through the product of the internal resistance R of the battery and the squared electrical current I:

q J = R • I 2 (II.10)
Deeper discussions are given on the choices for the parameterization of the internal resistance in the chapters describing the numerical solvers (see chapter III and applications in chapter IV).

II.3 Immersion cooling of Li-ion batteries II.3.1 Li-ion batteries in electric vehicles

Li-ion cells are gathered in modules to provide enough power to propel electric vehicles.

Each module is in a bunch of connected batteries. Usually, a module contains hundreds of cells. Then, several modules are connected and compose the whole battery pack powering the electric engine. The cells are connected following a pattern of series and parallel connections.

When designing a battery pack, the issue is to ensure a good energy density. The goal is to put the batteries close enough to each other to keep the overall volume as small as possible while respecting a minimal spacing, essential for thermal reasons explained thereafter. From a general perspective, the widespread use of electric vehicles for daily consumers impose specific constraints on the cars and battery packs manufacturers. The two main industrial objectives to meet are:

1. Provide high power engines capable of propelling normal size vehicles at high speeds.

2. Enable the fast charging of the battery pack. Encouraging the use of electric vehicles requires ultra-fast charging of batteries to close the gap between the time needed for filling the tank of a petrol car and charging an electric vehicle [START_REF] Tomaszewska | Lithium-ion battery fast charging: A review[END_REF].

These two capabilities induce a high electric current submitted to the packs. As seen earlier in the manuscript, these electric currents will provoke considerable heat loads on the batteries. Then, it is necessary to develop systems that allow the batteries to manage and retrieve the heat produced. The response to these constraints is the development of Battery Thermal Management Systems (BTMS).

During a classical operation of an electric vehicle, which consists of an alternation of charge and discharge sequences, the main thermal objectives to meet with the BTMS are: 2. Ensure a good temperature homogeneity between all the cells. The temperature of each cell should be as uniform as possible across the whole pack.

3. Inside a cell, keep a good temperature homogeneity. The temperature field inside a single cell should not present too significant variations and gradients.

Meeting those objectives is the solution to prevent the problems related to ageing and safety caused by temperature variations. Objective 1 tries to prevent the risk of the thermal runaway (maximal temperature) and not degrade the batteries' performance (minimal temperature).

For a given electric vehicle, it is preferable that the batteries age at the same rate and then present characteristics as homogeneous as possible. This fact is the primary motivation of objectives 2 and 3.

II.3.2 Battery thermal management systems (BTMS)

This BTMS usually includes a cooling element (liquid, gas, phase change material), conduct or enclosure to store the cooling element, a pump, a radiator and an expansion vessel (for fluid systems). BTMS are fully integrated into the vehicles, so the weight and volume of the whole system have to be considered in order not to penalize the power/weight ratio of the car.

Let's review more precisely the principle of the four leading technologies [START_REF] Al-Zareer | A review of novel thermal management systems for batteries[END_REF] available in the market. Phase change materials A phase change material (PCM) can present a phase change between solid/liquid within a narrow range of temperatures. The principle of BTMS using PCM is to insert the batteries within the specific material, which stores the heat produced by the batteries. It is seen as a passive technique, as no external energy has to be provided to the cooling medium to perform the heat transfer.

Cooling Circuit Battery Module

The batteries first generate the heat, the PCM retrieves this energy thanks to the phase change. Then the heat is diffused to the battery case by convection (natural or forced), in contact with ambient air cooling channels [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF]. The heat transfer process is illustrated in The advantage provided by the phase change is that the energy transfer operates at a quasi-steady temperature, the latent phase change temperature. Then this technique helps to get a good temperature uniformity between the cells in the pack [START_REF] Lazrak | An innovative practical battery thermal management system based on phase change materials: Numerical and experimental investigations[END_REF]. Practically, the material used is often paraffin. It presents adequate properties such as a narrow melting temperature range (between 30 and 60 o C ), a high latent heat, small volume variation during phase transition, stability, non-flammable and non-explosive [START_REF] Hallaj | A Novel Thermal Management System for Electric Vehicle Batteries Using Phase-Change Material[END_REF][START_REF] Rao | A review of power battery thermal energy management[END_REF].

c
Figure II.13: PCM cooling system with air channels and heat transfer principle at the cell level. Figure adapted from [START_REF] Al-Hallaj | Safety and Thermal management for li-ion batteries in transportation applications[END_REF] and [START_REF] Lazrak | An innovative practical battery thermal management system based on phase change materials: Numerical and experimental investigations[END_REF].

The following techniques are seen as active methods, involving the flow of a cooling fluid (liquid or air) within the battery pack. The flow is usually generated by a pump providing the energy necessary to ensure the forced convective heat transfer.

Indirect liquid cooling

In the indirect cooling BTMS, the cooling fluid is not in direct contact with the pack's batteries. Usually, a small channel meandering through the whole battery pack is used to ensure the flow of the coolant. Several configurations exist [START_REF] Ahmad A Pesaran | Battery Thermal Management in EVs and HEVs: Issues and Solutions[END_REF]:

tubes can be directly attached to cells or cooling plates can be disposed around the modules.

The plates are usually composed of several channels with a coolant flowing in each channel.

In The route of heat exchange in such configurations is the following: the heat is generated in the batteries. Then, the heat is spread through the channel material through conduction.

The liquid flowing in the channels is retrieving the energy by forced convective heat transfer.

One of the advantages of indirect cooling is its practicability to implement in electric vehicles configurations. The system presented by [START_REF] Cao | Thermal modeling of full-size-scale cylindrical battery pack cooled by channeled liquid flow[END_REF] in Fig. II.14c has a quite complex shape, but it integrates well in the overall car design. The necessary amount of coolant is not that large to fill the total volume of the channels. It is interesting regarding the weight of the whole system and thus the car's weight. The indirect cooling approach presents some drawbacks, especially for the temperature uniformity criteria. The circuit of the cooling channels across the pack is a determinant on the thermal performance of the system. The computations performed by [START_REF] Dubey | Direct Comparison of Immersion and Cold-Plate Based Cooling for Automotive Li-Ion Battery Modules[END_REF] and illustrated in Immersion cooling systems can be used with air or dielectric liquids as fluid performing the heat transfer. The cooling performances of systems using liquids are better than air [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF][START_REF] Ahmad A Pesaran | An Approach for Designing Thermal Management Systems for Electric and Hybrid Vehicle Battery Packs[END_REF]. However, using air allows to reduce the technical issues: the air can be exhausted directly out of the vehicle and replaced by fresh air using existing technologies in conventional thermal cars.

Immersion cooling systems are also promising to ensure a good temperature homogeneity in the modules and each cell. Still, the flow path of the coolant should also be carefully designed. Practically, on the one hand, the goal is to maximize the contact area between the fluid and the surface of the cells. Thus the heat transfer between the fluid and solid is maximized. On the other hand, the time spent between the inlet and outlet of the module should be reduced for a given particle of fluid. The fluid should meander between the pack cells to retrieve as much heat as possible. But when its temperature reaches a specific value, the fluid should be exhausted and replaced by 'new' and cold fluid.

To fulfill these two constraints, two parameters can be considered for a given fluid. First, the flow rate of the fluid crossing the module could be increased significantly. This would be efficient, but in practice, the energy consumed by the external pump can deprecate the energy balance of the whole system. Focus can also be set on the geometry of the module, fluid intake and exhaust, and the cells arrangement. A setup like the one in Fig. II.17 from [START_REF] Huat | Computational fluid dynamic and thermal analysis of Lithiumion battery pack with air cooling[END_REF] is interesting in that sense. The cooling air is entering the module above the aligned battery cells (Fig. II.17a). Then it travels the module from the bottom to the top before being exhausted in the superior conduct. This geometry and flow route ensures an excellent thermal homogeneity: i) in the whole pack, ii) within each cell. Indeed, for condition i), each cell receives inlet air at the same temperature as the flow is entering by the bottom uniformly.

For condition ii), the air is performing the convective heat transfer only along with the height of one cell. This short path for the traveling air does not spend too much time around the warm cells, and an even cooling has more chance to be achieved.

(a) Geometry of the immersion cooling BTMS. The immersion cooling systems are also efficient to prevent thermal runaway. An experiment has been conducted at Exoes [START_REF] Rémi | Novel battery thermal management enabling near zero temperature gradient for fast charging while improving safety[END_REF] for a pack of eight cells closely spaced, immersed in a high heat transfer coefficient dielectric fluid. One of the cells has been mechanically abused to initiate the thermal runaway. The experiment showed that the damaged temperature has considerably increased, as expected. However, the presence of the fluid between the cells prevented the temperature from spreading towards the other cells. Thus thermal runaway propagation in the battery pack has been avoided.

The drawbacks of immersion cooling systems are mostly related to their practical implementation in vehicles. They induce heavy machinery, and the design of a pack capable of handling the coolant flow with its pressure variations is technically challenging to address. We recall that this configuration is under development at Exoes, which we simulate with the two solvers presented later in this thesis. To summarize the lastly reviewed thermal management systems, the Tab. II. [START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF], [START_REF] Al-Zareer | A review of novel thermal management systems for batteries[END_REF] and [START_REF] Rémi | Novel battery thermal management enabling near zero temperature gradient for fast charging while improving safety[END_REF].

Exoes company aims to develop battery packs as the ones shown in

Two phased liquid immersion cooling

II.3.3 Physics of immersion cooling

To address the simulation of immersion cooling systems, the physical problem to solve is the case of conjugate convective heat transfer. This section will cover the physics and equations governing the conjugate heat transfer process between a fluid and a solid. Both forced and natural convective heat transfer are of interest in the study of immersion cooling systems. Forced convection is known to show increased cooling performances for Li-ion batteries thermal management perspectives. Also, the natural convection is studied for scenarios of thermal runaway, to assess the cooling capabilities of the system when the pump does not generate the flow rate anymore.

II.3.3.1 Conjugate heat transfer

The conjugate heat transfer consists of the coupling of the heat equation in a solid with the energy and mass conservation equations in the fluid domain surrounding the solid [START_REF] Horvat | Calculation of conjugate heat transfer problem with volumetric heat generation using the Galerkin method[END_REF]. It is used to compute the convective heat transfer between the heated solid (volumetric heat source term) and the fluid flowing around. This approach can be seen as an extension of the convective heat transfer computation using Newton's law of cooling [START_REF] Çengel | Heat and mass transfer: A practical approach[END_REF].

Convective heat transfer with law of cooling

The convective heat transfer can be represented using this phenomenological law. Then the heat transfer qconv in [W • m -2 ], is computed using:

qconv = h(T s -T ∞ ) (II.11)
where h is the convective heat transfer coefficient, T s the temperature on the solid surface, However, this first approach can be seen as a low fidelity way to compute the heat transfer between the solid and the fluid. It uses a phenomenological law of cooling and omits spatial dependencies in the fluid and solid domain. The first solver presented in this thesis (in chapter III), which is developed at Exoes, offers a similar vision to compute the heat exchanged between the solid and the fluid.

T

Conjugate heat transfer with CFD

The conjugate convective heat transfer approach has become practically possible with the rise of computational methods and resources [START_REF] Perelman | On conjugated problems of heat transfer[END_REF]. This model allows substituting the latter linear relationship between the temperature and the heat transfer rate. The temperature field of both fluid and solid domains and the associated heat transfer can be computed by skipping the evaluation of the heat transfer coefficient h. Each medium is seen as an independent domain, where the physical processes occurring in fluid, and solid subdomains are computed in each domain respectively, and coupled through continuity conditions. The heat transfer between the two media is reproduced by the conjugate conditions applied at the interface of the two objects. It comes to enrich the approach described by Newton's law of cooling as the temperature fields can be explicitly computed. The convection heat transfer rate can also be assessed at every point of the domain.

The conjugate heat transfer problem can be fully represented through partial differential equations in the fluid and solid domains. Let's write down the equations used to solve the conjugate heat transfer, in the case of a laminar and incompressible fluid flowing around a solid submitted to a volumetric heat source (Fig. II.20). To stick with the applications covered in this thesis, we consider that the fluid does not receive heat from a volumetric heat source. Note that further explanations are given in chapter IV on the physical assumptions and depending expression of the equations. However, we write already the equations here to facilitate the explanations on the physics at stake.

The unknowns are the temperature in the fluid and solid domains, the fluid velocity and pressure, denoted in Tab. II.5. The position in the domains is given by the coordinates stored in the vector x and the time is denoted t.

T s ≡ T s (x, t) Temperature field in solid domain

T f ≡ T f (x, t) Temperature field in fluid domain u ≡ u(x, t) Velocity field of the fluid p ≡ p(x, t)
Pressure field of the fluid Table II.5: Notations of the unknowns in the conjugate heat transfer problem.

In the fluid domain Ω f , the equations to solve are given in Eq. (II.12). The mass conservation is ensured through the incompressible Navier-Stokes equations. The temperature field of the fluid T f (x, t) is computed using the heat equation. Here ρ f is the density of the fluid, C p,f the specific heat, ν f the kinematic viscosity of the fluid and λ f the vector of thermal conductivities. F v represents the external forces applied to the fluid.

               ∇ • u = 0 ∂u ∂t + (u • ∇)u = ∇ • (ν f ∇u) - 1 ρ f ∇p + F v ρ f C p,f ∂T f ∂t + u∇T f = ∇ • (λ f ∇T f ) (II.12)
In the solid domain Ω s the temperature field T s (x, t) is computed through the unsteady conduction equation Eq. (II.13). ρ s is the density of the solid, C p,s its specific heat and λ s the vector of thermal conductivities. The volumetric source term heating the solid is denoted q g .

ρ s C p,s ∂T s ∂t = ∇ • (λ s ∇T s ) + q g (II.13)
At the interface, denoted Γ i , the conjugate coupling conditions traduce the temperature and heat flux continuity between the two domains. The conditions are expressed in Eq.

(II.14), involving the temperature evaluated on the interface, the vector normal to the inter-face n and the thermal conductivities in each medium.

T f,i = T s,i λ f ∂T f (x, t) ∂n = -λ s ∂T s (x, t) ∂n (II.14)
We illustrate schematically the process in Fig. II.20. Physically, the heat is produced in the solid domain through the source term q g . The heat is spread through the solid by the thermal conductivity λ s . Then, the heat transfer occurs between the fluid and solid thanks to the coupling conditions: the fluid layer around the solid edge retrieves the produced heat.

The behavior of the velocity field depends on the convective heat transfer regime considered:

forced or natural convection.

Forced convection

In the case of forced convection, a velocity is imposed at some boundaries of the fluid domain. The fluid flows around the heated object by respecting the no-slip condition at the interface. This method is very effective in terms of cooling performance.

Most of the immersion cooling systems mentioned earlier in this chapter use this process.

In a forced convection regime, the fluid layer next to the solid walls is heated by the solid through conduction and coupling conditions. Thanks to the imposed velocity, this layer is regularly evacuated. A new one replaces the 'old' heated layer with a lower temperature (usually imposed at inlet boundaries). Then heat transfer is improved as the solid is constantly facing cold temperatures [START_REF] Çengel | Heat and mass transfer: A practical approach[END_REF].

Natural convection

In this case, no velocity is externally imposed on the fluid. The heated solid creates temperature gradients in the fluid domain following the conjugate heat transfer process. Those temperature variations cause density gradients in the fluid domain. Those gradients generate hydrostatic disequilibrium and create a flow around the solid. According to Boussinesq's approximation, the vertical components of the flow are predominant in this case. However, the fluid moves at a significantly slower speed, or with a lower flow rate, around the solid compared to forced convection. So the 'hot' layer around the solid might not be evacuated as rapidly, and the temperature difference between the two domains is decreased, which deteriorates the heat transfer rate.

II.3.3.2 Heat transfer and flow characterization

As mentioned just above, the performances of the heat transfer in a conjugate heat transfer configuration will depend on parameters coming from physical processes occurring in the fluid and the solid domains. Regarding the forced convective heat transfer, the following dimensionless numbers are helpful to have insights into the predominant physical phenomena.

The following explanations are given in the case where the solid produces internal heat and the fluid around is cooler. It is given by ν = ρ η with ρ the density and η the viscosity of the fluid.

Prandtl number Considering the conjugate heat transfer of a flow around a heated object, we can expect variations in the heat transfer rate depending on the flow regime as detailed in [START_REF] Çengel | Heat and mass transfer: A practical approach[END_REF]. For a laminar regime, the fluid layers in contact with the solid will retrieve the heat produced and diffuse it by conduction through the so-called thermal boundary layer. The thermal boundary layer is the region over the obstacle's surface where the temperature gradient normal to the surface is significant. The convective heat transfer coefficient is higher in a turbulent regime as the fluid particles are more 'mixed' than in a laminar regime. Fluid particles with significant temperature differences are more likely to encounter and then exchange heat, which increases the overall heat transfer coefficient.

"The random eddy motion of a group of particles resembles the random motion of molecules in a gas -colliding with each other after traveling a certain distance and exchanging momentum and heat in the process.", Yunus A. Çengel in [START_REF] Çengel | Heat and mass transfer: A practical approach[END_REF].

Depending on the flow regime, these energy dissipation processes are characterized by the Prandtl number P r. This dimensionless number gives the ratio of the molecular diffusivity of momentum and the molecular diffusivity of heat. It allows comparing the temperature profile with the velocity profile in the thermal and velocity boundary layers. For instance, in oils that are usually really viscous fluids, the heat diffuses more slowly than momentum. They are efficient to transfer energy through convection compared to conduction. The heat diffuses much more quickly than the momentum in liquid metals, as these materials usually present a high thermal conductivity. The energy transfer occurs mostly due to the conduction within the thermal boundary layer. In this case, the temperature profile is almost independent of the velocity profile as the heat diffuses much more quickly as the fluid flows.

Nusselt number Back to heat transfer perspective only, the Prandtl number is helpful to assess the Nusselt number N u value. In forced convective heat transfer, this dimensionless quantity compares the predominant heat transfer between the fluid and a solid layer: convection or conduction. Locally at a point on the surface, it is expressed as:

N u = hL λ (II.17)
h is the coefficient of heat transfer. L is the characteristic length of the object, as used for the Re number computation. λ is the thermal conductivity of the fluid.

If the Nusselt number has great values (N u 1), the heat transfer by convection is more important than the conduction. The heat exchange is mostly due to the flow of the cold fluid retrieving the heat produced by the solid. When the Nusselt number presents small values (N u 1), the heat transfer between the solid and the fluid is mainly due to conduction. For instance, in the case of natural convection, the velocity of the fluid around the obstacle is low.

The heat is retrieved in the thermal boundary and dissipated through the fluid with thermal conductivity.

Practically, the Nusselt number can be used to compute the actual heat transfer coefficient of the fluid. In forced convection, for steady flow and fluid with constant thermal and mechanical properties, it can be linked to the Prandtl and Reynolds numbers with the Hilpert correlation [START_REF] Manohar | A Comparison of Correlations for Heat Transfer from Inclined Pipes[END_REF][START_REF] Frank | Fundamentals of heat and mass transfer[END_REF]:

N u = cRe α P r β (II.18)
The coefficients c, α, β are determined empirically depending on the specific case considered.

The common objective in thermal and fluid dynamics study is to assess the value of the Nusselt number through the computation of Re and P r, to finally get the value of the global (integrated on the surface of the object) or local value of the heat transfer coefficient h.

In natural convection, the Nusselt number can be assessed using the Rayleigh number Ra, computed through the Grashof number Gr under the Boussinesq' approximation. The Grashof number characterizes the intensity of the free convection occurring around the heated solid.

N u = cRa γ (II.19)
The choice of the coefficients c and γ depends on the type of flow generated by the convection, laminar or turbulent [START_REF] Frank | Fundamentals of heat and mass transfer[END_REF]. The Rayleigh number is given by:

Ra = P r • Gr
The Grashof number is computed with:

Gr = gβ T L 3 ν 2
The Grashof number is expressed through the temperature difference T along the surface of a vertical solid of height L, the acceleration of gravity g, the kinematic viscosity ν and the thermal expansion coefficient of the fluid β.

II.4 Chapter conclusion

This chapter has covered the main elements composing the Li-ion batteries. The physical parameters involved in the electrical characterization of batteries were detailed. The thermal issues related to batteries operation and the associated hazards were explained with physical considerations. The link between the heat generation and the electrical behaviour of the batteries was stated at this point. More details on this physical phenomenon were given by covering the modeling approaches commonly used in literature to represent the electrical processes in the batteries. The different scales considered could show the advantages and drawbacks of each strategy.

We also highlighted the need for Battery Thermal Management Systems for electric vehicles operation by detailing the thermal issues. The different systems existing in the industry were covered with their advantages and drawbacks. Specific attention was paid to the im-mersion cooling systems, being the physical case studied in the scope of this thesis.

Finally, narrowing the aim to the outcomes of the present work, we focused on the heat transfer modeling strategy selected in this thesis. First, the phenomenological approach, used in the low fidelity solver, was explained with the law of cooling. Then, regarding the high fidelity solver, the conjugate heat transfer approach was described by presenting the involved equations along with relevant physical considerations and parameters.

Overall, this chapter drew the landscape to motivate the scientific issues addressed in this work: develop accurate numerical models to represent the heat transfer of Li-ion batteries under immersion cooling heat transfer and take into account the uncertainty on physical parameters involved in the electrical and thermal behaviors of the batteries.

LOW FIDELITY SOLVER FOR IMMERSION COOLING: ICEXO

Chapter abstract

This chapter presents the first model, denominated ICExo. We first describe the equations solved by the thermal and electrical model. Then the model's performances are illustrated by comparison with an experimental test case. Using this test case, we illustrate the ability of the model to reproduce the measurements of temperature. Finally, some concluding remarks cover the advantages and the limits of the approach proposed with this numerical model. The heat transfer part is based on the resolution of the heat equation using the finite difference in a cylindrical battery seen in a 2D geometry. The thermal exchange between the fluid and the solid is computed, solving the equations of enthalpy balance through a conduct.

Outline

This approach allows considerable computational cost as the flow equations are not solved.

The model's accuracy is still ensured, as shown later in the chapter.

In

III.2 Thermo-electrical solver: ICExo

The model presented in this paper computes the transient heat transfer between a Li-ion battery, the surrounding fluid and the evolution of electrical parameters for the case of a battery submitted to unsteady input electric currents. The approach chosen was to develop a model allowing insights on some parameters of interest evolution in immersion cooling

III.2.1 Thermal equations

In this section, we describe the equations representing the heat transfer between the coolant

ρC p ∂T (x s , t) ∂t = λ∆T (x s , t) + q g , (III.1)
where ρ is the density of the Li-ion cell, C p its specific heat capacity and λ = (λ r , λ z ) is the vector of thermal conductivities in radial and axial directions respectively. q g is the volumetric source term representing the heat generated within the Li-ion cell by the Joule effect and is computed as follows

q ∝ R 0 • I 2 (t) (III.2)
where R 0 stands for the internal resistance of the cell, and I(t) is the electric current submitted to the batteries at a given time t. The heat is generated in the battery through the volumetric source term in Eq. III.2 and then is spread in the domain through Eq. III.1. As the computational domain represents only half of the Li-ion cell, a symmetry condition is applied for the temperature at the center of the battery Γ sym . At the interface Γ int between the fluid and solid domains, a Neumann condition is applied, giving the heat flux generated by the temperature difference between these domains:

∂T (x Γ int , t) ∂n = -h f S Γ int T (x f , t) -T (x s , t) (III.3)
where S Γ int is the surface of the whole boundary, T (x f , t) the temperature in the fluid domain and T (x s , t) the temperature in the solid domain at the other side of the interface. In practice, the solid domain is meshed, and a finite differences scheme is used for space and time to solve the heat equation Eq. (III.1). The mesh containing 100 cells is visible in Fig. III.2. It is meshed using a cartesian grid with the same number of cells in both radial and vertical directions.

The mesh is then more refined in the radial direction.

ṁ h(x f , t) -h(x f , t) = ϕ th (III.4)
where ṁ is the mass flow rate of the fluid in the conduct. The heat flux ϕ th defined in following represents the heat exchanged between the battery and the surrounding fluid in the control volume V :

ϕ th = h f • S V T (x f , t) -T (x s , t) (III.5)
where h f is heat transfer coefficient of the fluid. S V is the surface of the interface between the solid domain and the control volume.

At the fluid inlet boundary Γ inlet , a Dirichlet condition is applied for the specific enthalpy:

h(x Γ inlet , t) = h ini . (III.6)
At the fluid outlet boundary Γ outlet , an homogeneous Neumann condition is applied:

∂h(x Γ outlet , t) ∂n = 0 (III.7)
The temperature of the fluid is obtained by interpolation of tables linking the specific enthalpy, the pressure, and the temperature using the REFPROP database [START_REF] Lemmon | NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 10[END_REF]. In the following, the interpolation operator of 2D gridded data for an output quantity u is written as:

u = Interp 2 (x q , y q , [T1, T2, TU])
With this notation, the table [T1, T2, TU] contains the values of the quantity to compute u ∈ TU, depending on the values in the vectors T1 and T2. The value u is computed for the query points x q ∈ T1 and y q ∈ T2, by interpolation on the values in TU corresponding to the query points. Numerically speaking, the Interp 2 function utilized is the Matlab 'interp2' routine (table lookup). The value u is obtained via linear interpolation of the data available at the neighborhood of the query point (x q , y q ).

The temperature of the fluid is then linked with the pressure and enthalpy in the table

written [h, P, Temp]. Then, temperature of the fluid at a location x f and a time t is given as follows: The specific enthalpy h ini at the inlet fluid boundary is given by the following interpolation:

T (x f , t) = Interp 2 h(x f ,
T ini = Interp 2 h ini , P 0 ; [h, P, Temp] REFPROP (III.9)
where T ini is the temperature at the beginning of the simulation t 0 and also stands for the steady temperature outside the computational domains. . This temperature is chosen to simulate the battery's skin temperature as if it was located within a gap between two cells in a more extensive pack configuration.

This choice has been encouraged by comparing the numerical results with the experimental data, shown later in this chapter. To get a scalar quantity, the mean of the temperatures in this area is computed, giving at any time t the temperature denoted T qoi (t). At any time t, the temperature of interest T qoi (t) is computed using

T qoi (t) = 1 S Ω right int Ω right int T (x, t)dx.
(III.10)

III.2.2 Electrical equations

Along with the thermal equations, electrical parameters of the Li-ion cells are also computed during the simulated time. This study focuses on the computation of the voltage and state of charge evolution of the Li-ion battery.

The state of charge SOC of the cell is computed at each time t of the simulation using the relation:

SOC(t) = 1 • SOH - I(t) C • 3600 • t (III.11)
where C stands for the capacity of the cell in to the voltage source of the battery E m and its internal resistance R 0 . Following this approximation, the full voltage V of the battery is given by the following relation:

[A • h]. SOH
V (t) = E m (t) + Z(t) • I(t). (III.13)
The voltage source of the circuit is the electromotive force denoted E m . Z is the impedance of the RC circuit representing the battery. E m represents the no-load voltage of the open circuit. It is computed at each time t of the simulation using the following relation:

E m (t) = E i (t) -E T (t). (III.14)
Note that E i is the voltage of a cell with no input electric current, and its value is obtained by interpolation of tables linking the electrical current, the state of charge and the voltage, as follows

E i (t) = Interp 2 I = 0, SOC(t); [I, SOC, E i ] . (III.15)
Note also that E T is the voltage of the cell dependent on the initial temperature and is obtained by interpolation of tables linking the temperature, the state of charge and the voltage, computed as follows

E T (t) = Interp 2 T ini , SOC(t); [Temp, SOC, E T ] . (III.16)
The table E i and E T contain values of the voltage with average discrepancy of 0.

Z(t) = R 1 (t) 1 + R 1 (t)C 1 (t) + R 0 , (III.17)
where R 1 and R 0 are the internal resistance of the RC circuit. R 0 represents the main contribution to the resistance of the battery. A constant value is chosen for this parameter.

The resistance R 1 is obtained by interpolation of tables linking the temperature, the state of charge and the resistance: Note that in practice, R 0 is assumed much more significant than R 1 . To simplify the thermal model behavior and in the perspective of performing the calibration of the internal resistance parameter, the choice was made to consider only the contribution from R 0 for the internal resistance in the heat generation of the Li-ion battery in Eq. (III.1). The battery manufacturer provides all the data tables used in the previous equations. In addition, this model allows simulating the electrical parameters for an equivalent battery module rather than a single cell. In this context, the type of connections in the battery pack is considered, i.e. the parallel and serial configuration of the connections between all

R 1 (t) = Interp 2 T ini ,
Li-ion cells. The equivalent input electrical current and voltage can be computed using Kirchoff's laws. For the thermal part, the values of the parameters involved in the heat transfer, like the flow rate, are usually given for the whole module. To go from the vision of the whole pack to one equivalent cell, these parameters are divided by the number of cells in the module considered.

III.3 Validation of the model with an experimental test case

In this section, an original experimental test case conducted by Exoes is presented. It aims to study the thermal behavior of a small pack of batteries immersed in a coolant. The batteries are heated under various charging/discharging conditions. The purpose of the study is to produce trustful and reliable experimental measures of temperature, voltage and current for this small immersion cooling configuration. This kind of data interests anyone trying to calibrate numerical models for such problems. The experimental results are provided in a dataset available on a Mendeley Data repository [START_REF] Solai | Lithium-ion Batteries Immersion Cooling Experiment[END_REF].

Once the full details on the experimentation are given, we illustrate the accuracy of the model ICExo with the reproduction of the experimental test case.

III.3.1 A review of experimental setups for immersion cooling

Many experiments of immersion cooling configurations have been conducted in the literature.

While experimental facilities and methods can differ, measurements are often related to the temperature evolution of Li-ion cells under various conditions. Notably, most references study the effect of charge and discharge rates on the batteries temperature evolution. These facilities usually feature a pack composed of few cells, going from a lonely cell as in [START_REF] Sundin | Thermal Management of Li-Ion Batteries With Single-Phase Liquid Immersion Cooling[END_REF] up to larger configurations as in [START_REF] Fan | Experimental study on the thermal management performance of air cooling for high energy density cylindrical lithium-ion batteries[END_REF] The immersion cooling experimental setups usually allow studying the impact of the geometry and the arrangement of the cells. The references [START_REF] Tan | Numerical investigation of the direct liquid cooling of a fastcharging lithium-ion battery pack in hydrofluoroether[END_REF][START_REF] Wang | Calculation methods of heat produced by a lithium-ion battery under charging-discharging condition[END_REF][START_REF] Fan | Experimental study on the thermal management performance of air cooling for high energy density cylindrical lithium-ion batteries[END_REF] studied the influence of arrangement geometries on the cooling effectiveness. Indeed, [START_REF] Fan | Experimental study on the thermal management performance of air cooling for high energy density cylindrical lithium-ion batteries[END_REF] showed experimentally that an aligned arrangement in straight rows of cells displays the best cooling capabilities, based on temperature measurements within the pack. The thermal regimes cooling performances are also studied through those experimental setups. [START_REF] Li | Experimental Investigation on a Thermoelectric Cooler for Thermal Management of a Lithium-Ion Battery Module[END_REF] assess the performances of a newly developed cooling system, a thermoelectric cooler, by comparing it against natural and forced convection cooling using the same facility. [START_REF] Trimbake | Mineral Oil Immersion Cooling of Lithium-Ion Batteries: An Experimental Investigation[END_REF] also assess capabilities of their innovative oil impingement cooling technique. Some authors use experimental setups to build or calibrate numerical models. For instance, [START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF] assess the heat generated by the batteries from temperature measurements and include it as a source term in the model to ensure good numerical reproducibility of the experimental data. Many cited experiments use steady discharge rates in each investigated experimental condition. Only a few references perform temperature measurements while alternating charge, discharge, and rest cycles at different rates. In [START_REF] Sundin | Thermal Management of Li-Ion Batteries With Single-Phase Liquid Immersion Cooling[END_REF][START_REF] Trimbake | Mineral Oil Immersion Cooling of Lithium-Ion Batteries: An Experimental Investigation[END_REF] the temperature evolution the cells are investigated during such cycles.

The present experiment gathers many of the above features, providing novel experimental data. Precisely, the setup consists of a lab-scale battery pack of eight 18650 cylindrical cells immersed in a coolant and under a forced convection regime. The spacing between each cell and with the enclosure walls is narrow. In that sense, the current setup gets closer to realistic geometries of manufactured battery packs. Thermocouples measure temperatures at cathode and anodes locations for some of the eight cells. The experiment aims to measure the temperature evolution of the batteries submitted to different charging and discharging cycles.

In addition, electrical characteristics such as voltage and electrical current are monitored and put in perspective with the temperature evolution. The several cycles investigated present various charging rates within the same experiment. The experimental measurements allow quantifying the effect of charging rates in a transient regime and how the cooling capabilities are affected by such variations.

III.3.2 Experimental test case: lab-scale immersion cooling

The experimental case setup is represented in Thus, the values obtained in the plots for the electrical current are eight times higher than the values specified in Tab. III.2 (given for one cell), according to the wiring of the pack.

For both cycles, the temperature discrepancy between the cathode and the anode temperatures is explained by their respective enclosure locations. The cathode of Li-ion cells is at the entrance of cooling fluid in the enclosure, while the anode is at the exit. Consequently, the cooling fluid retrieves some heat while flowing around the cells. Also, from t ≈ 50 [sec] in RC cycle and the beginning of the DS cycle, the fluid temperature at the inlet is lower than the temperature at the outlet, which agrees with the cooling process taking place within the enclosure.

Overall, the flow rate measurements show that the steady condition is fulfilled, even if the rate presents some small noisy variations.

III.3.3 Numerical results

The values for the thermal conductivities and the specific heat capacity to initiate the calibration procedure was inspired from the references [START_REF] Wei | A comprehensive study on thermal conductivity of the lithium ion battery[END_REF][START_REF] Drake | Measurement of anisotropic thermophysical properties of cylindrical Li-ion cells[END_REF].

Parameter Symbol Value Units

Mass flow rate ṁ 2.925

• 10 -3 kg • s -1 Heat transfer coefficient h f 2.140 • 10 +2 W.m -2 .K -1 Internal resistance R 0 3.253 • 10 -3 Ω State Of

III.3.4 Discussion on the numerical results

As shown above, the numerical model can reproduce the physical process of heat transfer between the heated batteries and the fluid in a forced convection regime. The heating rate change with respect to the cooling is also well captured (see results of DS cycle). Finally, the way the model was developed allows the simulations at a low computational cost.

However, the modeling approach performed in this model is based on significant assumptions and simplifications. There are features and aspects of the real physics occurring that the model is, by nature, not able to reproduce. In particular, the coolant flow is represented only by the mass flow rate. There is no data given by the model allowing to evaluate the coolant's velocity field. Also, the geometry is not really represented. The heat transfer is computed with an 'averaged' cell, and the energy balances between the relevant quantities. The spatial and temporal dependency of the temperature is expected to be impacted strongly by these effects. The next chapter addresses the latter issues using a CFD model solving the full system of partial differential equations to represent the conjugate heat transfer.

III.4 Chapter conclusion

In this chapter, we presented a fast and accurate model to simulate the thermal behavior of heated batteries immersed in a cooling fluid. The electrical parameters of the battery are also computed using this model. The originality of the approach proposed in this model, performing energy balances instead of solving the flow equations, computing parameters through tables interpolation, allows to save considerable computational resources, compared to what would be done in a more classical CFD approach.

An original experimental setup of immersion cooling was described, and the experimental results were analyzed, trying to correlate the observed electrical and thermal phenomenon.

These data were used to validate the code behavior under different physical conditions. Overall, the code showed an excellent capability to reproduce the physical behavior in the present test case. However, the assumptions made to construct the model presented some limitations when assessing the code predictions. The comparison with the experimental data allowed some hypotheses on the discrepancy between the numerical and experimental results.

This final discussion highlighted the need to go with a higher fidelity solver, using CFD, to reproduce with better accuracy the thermal behavior of Li-ion batteries under immersion cooling.

COMPUTATIONAL FLUID DYNAMICS SOLVER: TRIOCFD

Chapter abstract This chapter focuses on the description of the CFD solver used in this thesis, the code TrioCFD. We state the hypothesis to address the simulation of immersion cooling of Li-ion batteries, the set of governing equations and the numerical methods. Two test cases are presented to assess the code capabilities and to study the thermal behavior of immersion cooling of batteries. The first test case is a comparison code to code on forced convective heat transfer simulation. The second test case deals with the simulation of the natural convection heat transfer and the validation with respect to experimental data available in the literature.

Outline

IV.1.2 CFD for immersion cooling: technical literature review

For a deeper understanding of the modeling approach performed with TrioCFD on immersion cooling systems, let us review the typical outcomes and interesting technical issues addressed in some literature references.

A 2D CFD conjugate heat transfer model is validated with experimental results for an original thermal management method using reciprocating airflow in [START_REF] Mahamud | Reciprocating air flow for Li-ion battery thermal management to improve temperature uniformity[END_REF]. The module they try to simulate is significant as it includes more than 20 Li-ion cells. In practice, only eight cells were represented with the CFD solver thanks to many symmetries in the geometry along flow channels. The heat generation is computed at the battery level, using the source term as described in section II.2.2.2. They also consider a model for the internal resistance, depending on temperature, playing a role in the Ohmic heating part of the source term. Besides, a lumped thermal model is also employed to represent this cooling facility. They compared By using simpler geometries, analytical models can be developed for immersion cooling configuration as it is done in [START_REF] Chalise | Conjugate Heat Transfer Analysis of Thermal Management of a Li-Ion Battery Pack[END_REF]. The results of their analytical model present a good comparison with the finite element simulations. The use of the analytical model allows to investigate with small computational effort the link between the temperature of the batteries and geometry arrangement. Also, they can assess the impact of the C-rate on the temperature in the domain. So they compare the numerical and analytical results under this perspective in steady-state cases.

Regarding the geometry of the arrangement of the cell, the authors from [START_REF] Wang | Thermal investigation of lithium-ion battery module with different cell arrangement structures and forced air-cooling strategies[END_REF] perform a comprehensive study on this characteristic using a 3D CFD solver. A flow of incoming air cools the batteries. The position of the fan generating the flow is also investigated. They represent a significant amount of cells (around twenty) immersed in the airflow and then highlight some typical heat transfer patterns depending on the geometries considered. The cooling performances were demonstrated to be intensely dependent on the airflow direction and path. A deterministic optimization problem was also conducted to obtain the best intercell distance in cooling efficiency.

Regarding the impact of the source term on the thermal response, some authors, such as [START_REF] Huat | Computational fluid dynamic and thermal analysis of Lithiumion battery pack with air cooling[END_REF] and [START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF] use the information provided by experimental data to adjust the CFD model and get the temperature response close to the measurements. More specifically, [START_REF] Huat | Computational fluid dynamic and thermal analysis of Lithiumion battery pack with air cooling[END_REF] tried different heating regimes and assessed the rate of temperature rise for the batteries.

As covered in the section II.2.2 from chapter II, the scale considered to represent the heating processes due to electrochemical reactions is of relative importance to the simulation accuracy and the involved models and resources. All the studies presented above consider a heat source term at the cell scale (see section II.2.2.2). The model presented in this thesis uses this approach too. Only a few studies present a more comprehensive vision, solving the partial differential equations at the electrochemical level, coupled with the heat generation in the cell and finally with the heat transfer with the surrounding fluid. Such an approach can be found in these work [START_REF] Suresh Patil | A novel dielectric fluid immersion cooling technology for Li-ion battery thermal management[END_REF][START_REF] Jithin | Numerical analysis of single-phase liquid immersion cooling for lithium-ion battery thermal management using different dielectric fluids[END_REF], using the so-called NTGK model.

Finally, note that the CFD solver used in the above studies for the thermal and fluid equations is Fluent for most cases.

Our approach with TrioCFD

In this chapter, we first recall the complete set of equations to solve the conjugate heat transfer problem. The equations from chapter II are detailed by considering all the physical hypotheses stated in the following test cases. Then, the employed numerical methods are described. The turbulence model choice is also detailed when it is integrated.

The code is compared with other numerical and experimental results using two different test cases. First, the solutions of TrioCFD are compared to Fluent on a case of a 3D simulation of immersion cooling of a pack composed of 9 Li-ion batteries [START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF]. We illustrate the code's mesh and time step convergence in this case. The purpose of the test case is to validate the overall behavior of the model, i.e. the choice of the heat generation, the resolution of flow and energy equations in all domains, on a configuration of forced convective heat transfer.

Then, the code is compared against an experimental test case of literature [START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF]. Here the attention is focused on the natural convective heat transfer, and thus the focus is set on the importance of the heat source model in the batteries. The results of 2D and 3D simulations on the same test case are compared with the experimental data. This comparison allows a

IV.2 Governing equations and numerical methods in TrioCFD

IV.2.1 Description of the physical problem

In this section, the equations written in the brief description of the Conjugate Heat Transfer approach in II.3.3.1 are further detailed. The equations solved in the solid and fluid domains are described considering the physical assumptions and simplifications. More precisely, the simulations addressed with the CFD code TrioCFD in this thesis consider the features and hypothesis listed in Tab. IV.1:

1 Transient conjugate heat transfer between a heated solid and a surrounding fluid.

2 Newtonian fluid.

3 Incompressible flow.

4 Homogeneous viscosity, density, specific heat, thermal conductivity.

5 Isotropic thermal conductivity unless specified otherwise.

6 Natural or forced convection.

7 Heat source from the solid domain and no heat source from the fluid domain. Table IV.1: Physical and simulation features of the CHT problem.

In the following, equations are expressed locally for a position x in the solid domain Ω s , in the fluid domain Ω f or at the interface between the media Γ, for a time t ∈ [t 0 , t f ] of the transient simulation.

IV.2.1.1 Solid domain: heat equation

In the solid domain denoted Ω s , the energy conservation is solved, in each location x ∈ R d (d = 2 or 3) and time t. In its more general form, the balance of internal energy e reads:

∂ρ s e ∂t + ∇ • (ρ s u) = -∇ • F -p∇ • u + τ : ∇u + q g (IV.1)
The left-hand side represents the variation of energy by volume unit. The right-hand side represents the heat flux by volume unit (first term). The two other terms stand for the power of the pressure and viscous forces. The last term q g stands for volumetric heat source coming from external phenomena. The domain is solid in this specific case, so the local velocity u is null. Similarly we have ∇ • (ρ s u) = 0.

The variation of the internal energy can be expressed in terms of temperature in the domain, denoted as T s ≡ T s (x, t) and considering the specific heat capacity at constant pressure, C p,s . Also, the heat flux in the domain is expressed with the local Fourier's law, considering the uniform and isotropic thermal conductivity (λ s ∈ R ): F = -λ s ∇T s The density ρ s is also homogeneous and constant. Finally, the heat equation in the solid domain reads :

ρ s C p,s ∂T s ∂t = λ s ∆T s + q g (IV.2)
Physically, this equation states that the internal variation of energy with time is equal to the heat flux dissipated through the domain added to an internal heat generation (q g > 0) or reduction (q g < 0). In the present work, we consider that the Li-ion batteries are producing heat, and the internal heat source term will be taken positive.

IV.2.1.2 Fluid domain: Navier-Stokes and energy conservation equations

The equations solved in the fluid domain give the velocity u and pressure p at a location x and a time t. u ≡ u(x, t) p ≡ p(x, t)

(IV.3)
The first equation is the local mass conservation in the fluid domain Ω f .

∂ρ f ∂t + ∇ • (ρ f u) = 0 (IV.4)
As stated in the hypothesis above, the flow is assumed incompressible in the present work so this equation reduces to:

∇ • u = 0 (IV.5)
Also, the conservation of momentum is ensured in the fluid domain. Under its local and conservative form it reads:

∂(ρ f u) ∂t + ∇ • (ρuu) = -∇p + ∇ • τ + F v (IV.6)
In this case, we consider a constant density ρ f and viscosity η f and the tensor for the rate of pure deformations is symmetric. Then the constraints tensor τ reduces to:

τ = 2η f 1 2 (∇u + ∇ T u) = 2η f ∇u
As the flow is considered incompressible (∇ • u = 0) we finally obtain:

∇ • τ = η f ∆u
Finally, the momentum equation to solve in the fluid domain is expressed as:

ρ f ∂u ∂t + (u • ∇)u = η f ∆u -∇p + F v (IV.7)
For the external forces, the gravity is considered in the 3D simulations and then we have

F v = -ρ f g
The gravity is expressed through g = -ge z with e z the basis vector in the vertical direction. In the case of 3D computations with natural convection, the Boussinesq approximation induces another term for the external volumetric forces. Following this approximation and considering incompressible flow, the density is dependent on the temperature of the fluid only. The buoyancy effects are accounted only due to the gravitational forces. The Boussinesq approximation applies if the variation of density |δρ| is really small compared to the real density: |δρ| ρ f . Then the external forces are expressed using the thermal dilatation coefficient of the fluid β T as:

F v = -ρ f gβ T (T f (x, t) -T 0 ) (IV.8)
T f being the local temperature of the fluid and T 0 a reference temperature to define. Otherwise, no external force is applied to the fluid.

Energy conservation in the fluid domain

∂ρ f e ∂t + ∇ • (ρ f u) = -∇ • F -p∇ • u + τ : ∇u + q f (IV.9)
The flow is incompressible (p∇ • u = 0), the heat generated by viscous effect are neglected (τ : ∇u = 0) and there is no internal heat generation in the fluid domain (q f = 0). Furthermore, the internal energy is expressed through the temperature of the fluid T f using the density ρ f and the heat capacity at constant pressure C p,f . In the same way, the heat flux within the fluid domain is expressed through Fourier's law: F = -λ f ∇T f Finally, the energy conservation in the fluid to solve reduces to: The averaged terms u, p represent the statistical mean of the velocities and pressure at each

ρ f C p,f ∂T f ∂t + u∇T f = ∇ • (λ f ∇T f ) (IV.
ρ f ∂u ∂t + (u • ∇)u = η f ∆u -∇p -ρ f u u + F v (IV.14)
This formulation induces an additional term: the Reynolds tensor ρ f u u representing the turbulent agitation. To express this tensor, a low Reynolds (low-Re) turbulence model is chosen based on the k -ε approach: the Jones and Launder model [START_REF] Jones | The prediction of laminarization with a two-equation model of turbulence[END_REF]. Following the Boussinesq's hypothesis on turbulence, the Reynolds tensor is expressed as:

u u = µ T ρ f (∇u + ∇ T u) + 2 3 kI (IV.15)
I is the identity matrix. This expression of the Reynolds tensor introduces the turbulent viscosity µ T . This viscosity is expressed by integrating the turbulent kinetic energy k and ε the dissipation rate of turbulent kinetic energy. Following the Jones and Launder model, the turbulent viscosity is given by:

µ T = c µ f f µ f ρ f k 2 ε (IV.16)
After including those two variables, the system is finally closed by solving additional transport equations for the averaged k and ε:

ρ f ∂k ∂t + u∇k = ∂ ∂y µ f + µ T σ k ∂k ∂y + µ T ∂u ∂y 2 -ρ f ε -2µ f   ∂k 1/2 ∂y   2
(IV.17)

ρ f ∂ε ∂t + u∇ε = ∂ ∂y µ f + µ T σ ε ∂ε ∂y + c 1 f 1 µ T ε k ∂u ∂y 2 (IV.18) -c 2 f 2 ρ f ε 2 k + 2µ f µ T ∂ 2 u ∂y 2
The values of the constant coefficients c µ f , c 1 , c 2 , σ ε , σ k and f 1 , f 2 , f µ f are given according to the original reference of the Jones-Launder model in [START_REF] Jones | The prediction of laminarization with a two-equation model of turbulence[END_REF].

IV.2.1.3 Coupling conditions at the interface

We assume the thermal contact is perfect at the interface between the solid and fluid domains.

So the temperature and heat flux are continuous. These coupling conditions ensure the heat transfer between the two media. So for each location on the interface x ∈ Γ and a time of the simulation t ∈ [t 0 , t f ], the coupling conditions are expressed with:

T f (x, t) = T s (x, t) Φ f (x, t) = Φ s (x, t) ⇐⇒ λ f ∂T f (x, t) ∂n = -λ s ∂T s (x, t) ∂n (IV.19)
The flux continuity is computed using the above formula, involving the solid and fluid thermal conductivities λ s , λ f and the vector normal to the interface n.

IV.2.2 Numerical methods

Let's detail the methods used for the spatial and time discretization to solve the transient conjugate heat transfer problem.

IV.2.2.1 Spatial and time discretization

Spatial discretization

The CFD simulations presented in the framework of this thesis are performed using an original Finite Element Method, introduced by [START_REF] Emonot | Méthodes de volumes éléments finis : applications aux équations de Navier-Stokes et résultats de convergence[END_REF]. This method uses elements derived from the Crouzeix-Raviart element [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations I[END_REF][START_REF] Heib | Nouvelles discrétisations non structurées pour des écoulement de fluides à incompressibilité renforcée[END_REF][START_REF] Thomas | Une méthode éléments finis à décomposition L2 d'ordre élevé motivée par la simulation de'écoulement diphasique bas Mach[END_REF]: nonconforming Crouzeix-Raviart elements. The method is employed on triangular meshes in 2D and tetrahedra in 3D.

The elements are described in [START_REF] Angeli | FVCA8 Benchmark for the Stokes and Navier-Stokes Equations with the TrioCFD Code-Benchmark Session[END_REF]: the pressure is discretized on the primary grid while the velocity is defined on a staggered dual grid. The degrees of freedom for the velocity are located at the center of the faces for each element. This element is convenient as it directly ensures the null divergence for the velocity. were not used here. We use the implicit solver based on the PISO method (Pressure Implicit with Split Operator) [START_REF] Guermond | On stability and convergence of projection methods based on pressure Poisson equation[END_REF]. The solver is running with the GMRES method and a Cholesky matrix preconditioner to solve the full linear system.

IV.2.2.2 Numerical coupling for conjugate heat transfer

Numerically, using the discretization described above, the coupling conditions ensuring the conjugate heat transfer are implemented by imposing the temperature equality on the center of the faces for the mesh cells at each side of the interface Γ. Considering this methodology, the mesh cells must be coincident all along the interface boundary.

Considering the time schemes, a parallel coupling strategy is employed [START_REF] Radenac | Methodology of numerical coupling for transient conjugate heat transfer[END_REF]. The strategy is summed up in

λ f ∂T f (x, t N ) ∂n Γ = -λ s ∂T s (x, t N ) ∂n Γ (IV.20)
Once the temperature at the interface is obtained, it is applied as a boundary condition on

IV.3 Simulation of 3D forced convective heat transfer

This section aims to show the capabilities of the code TrioCFD on an immersion cooling problem with forced convective heat transfer. This test case is reproduced from the reference [START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF]. More precisely, the results of this test case try to show the three following features. First, prove the numerical convergence of the code (mesh and time discretization). Then using the test case, we try to get insights into the physical phenomenon occurring in the immersion cooling configuration. Finally, we compare the results of TrioCFD with the present reference.

This code to code comparison highlights the specificities and the important features of this type of problem.

The problem covered here is a simplified representation of a Battery Thermal Management System. This study represents only a battery module composed of nine cells immersed in an airflow to make computational experiments affordable. Each cell generates heat due to the galvanostatic discharge current and its internal resistance. The forced convection regime then allow the airflow to retrieve the energy generated by the cells. Here, the battery module is represented as a part of a bigger battery pack containing more cells. The reader can refer to [START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF] to see more precisely the geometry reduction that has been made here.

First, we describe the geometry, the test case's physical properties, and the hypothesis made for the resolution. Then the results present first the mesh and time step convergence.

Then some plots of the velocity and temperature fields show the typical patterns of conjugate heat transfer problems. Finally, the results of TrioCFD are added to the plots of the current reference to give a proper code to code comparison.

IV.3.1 Test case description IV.3.1.1 Geometry

The battery module is composed of nine cylindrical cells. The radius of each cell is R = 21.2

[mm] and the height is L = 97.

IV.3.1.2 Physical parameters

An inlet velocity is imposed at the entrance of the fluid domain to generate the flow of air and then the forced convective heat transfer. Then we solve the conjugate heat transfer here as described in the previous section of this chapter. Note that here no turbulence is added to the Navier-Stokes equations to reproduce the approach performed by [START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF]. Based on the physical data described right after, the Reynolds number for this test case is around Re ≈ 2.7 • 10 2 .

The physical parameters of the batteries and the air are given in Tab. IV.2.

IV.3.1.3 Source term in the heat equation

The principle of the test case is to heat the Li-ion cells with a galvanostatic current during a finite time of 500 seconds. Each cell is submitted to a constant electric current of 23.36 [A] from t 0 = 0 up to t f = 500 sec. The simulation starts with an initial temperature in fluid and solid domains of

T ini = 295.1[K].
In [START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF]

Air properties

ρ f 1.2 kg.m -3 C p,f 1.0047 • 10 +3 J.kg -1 .K -1 λ f 0.0262 W.(mK) -1 η f 1.8426 • 10 -5 kg.(m.s) -1
Table IV.2: Geometric and physical properties of the battery cells generation source term q g is divided into two terms. It is precisely the approach described in the section II.2.2.2. One term represents the heat generation due to current transfer across the internal resistance induced by the electrochemical reactions in the cell. This term is the so-called Ohmic heating, denoted here q J . To this term is added the entropic heating, denoted q e . This heat source comes from the entropy changes caused by the electrochemical reactions during the simulation time.

Then the expression of heat source term for the heat equation Eq. (IV.2) reads:

q g = a (q J + q e ) = a      I(E oc -U ) ohmic heating -IT dE oc dT entropic heating      (IV.21)
a is a coefficient to match the physical dimensions between the source term and the heat equation (surface or volume normalization).

For the joule heating, a more applicable expression related to the known parameters of the Li-on cells is given by:

q j = R • I 2 (IV.22)
The internal resistance of the cell is expressed through the parameter R. In this study, an explicit dependence of R with the temperature T of the Li-ion cell is given by the relation:

R(T ) = -1.0 • 10 -4 • T 3 + 1.34 • 10 -2 • T 2 -5.345 • 10 -1 • T + 12.407 (IV.23)
where the temperature T = T s (x, t) is the temperature solved in the solid domain, expressed in Kelvins. In this relation the resistance is in [mΩ]. This formula was proposed by [START_REF] Park | Dynamic Thermal Model of Li-Ion Battery for Predictive Behavior in Hybrid and Fuel Cell Vehicles[END_REF] based on experimental measurements and has been used in several references.

IV.3.1.4 Boundary conditions

The boundary conditions applied to the solid and fluid domains are given in Tab. IV.

IV.3.2 Results

IV.3.2.1 Numerical convergence

Mesh convergence First of all, we want to ensure the grid convergence of the present computation. Six meshes composed of an unstructured grid of tetrahedra have been generated.

No specific refinement criteria have been applied at the solid and fluid domains interface.

This decision is made to stick with the reference settings used for the comparison. The mesh is generated uniformly regardless of the domain considered (fluid or solid). The simulation is unsteady, so the quantity of interest to assess the convergence is the temperature evolution of

IV.3.2.2 Physical analysis of the results

First of all, we visualize the transient temperature response of the cells. In More precisely, we can discuss the temperature homogeneity inside each cell. In Fig. IV.9, the temperature field in the solid domain is plotted together with the velocity vector field. About the airflow, the air velocity is very low behind each cell. There is a recirculating area behind the cells. In these zones, the cells present hot spots with higher temperatures. For instance, the temperature on the surface of cell 5 shows a difference of around 1 degree between the bottom and top of the cell. To get further insight on the link between the speed of air and the temperature hot spots, we show in 2D plots the temperature fields of solid domain and vector fields of air velocity, in Then following the Hilpert correlation at such flow regime (40 ≤ Re ≤ 4000) , the Nusselt number is computed with the following relation [START_REF] Manohar | A Comparison of Correlations for Heat Transfer from Inclined Pipes[END_REF]:

N u = 0.683 • Re 0.466 • P r 1/3
We obtain in this case N u = 8.4. Then, we can consider the heat transfer by convection is more significant than conduction in this configuration.

IV.3.2.3 Code to code comparison

To get a better insight into the quality of the results reproduction, the authors of [START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF] show the temperature values along some lines within the fluid domain ( About the direct comparison between the codes, the temperature rise in Fig. IV.12 is a bit overestimated by TrioCFD for a given height. However, this discrepancy between the codes is consistent with the velocity plots. FLUENT seems to underestimate the flow speed, which could explain the higher temperatures predicted.

IV.3.3 Discussion

This simulation highlights that this configuration can present some local hot spots due to a bad efflux of the coolant fluid. Specific attention should be paid to these zones when designing a BTMS. The cells arrangement should be designed to maximise the liquid velocity throughout the battery module. For instance, this objective can be achieved by spacing the cells as much as possible or by using a staggered geometry like it is explicitly studied in some references numerically [START_REF] Wang | Thermal investigation of lithium-ion battery module with different cell arrangement structures and forced air-cooling strategies[END_REF] or experimentally [START_REF] Fan | Experimental study on the thermal management performance of air cooling for high energy density cylindrical lithium-ion batteries[END_REF].

Both of the codes presented results with a good overall agreement. Some variations stand in the temperature predictions at some specific locations. The difference in the flow predictions can explain these variations in the temperature predictions. To gain confidence in the capabilities of TrioCFD to address this type of simulation, we present in the next section a comparison with experimental measurements of batteries temperature. This next case focuses on results by natural convection to simplify the thermal problem, reducing the velocities of air involved in the resulting thermal predictions.

IV.4 Simulation and validation of natural convection heat transfer

In this section, we validate the results of TrioCFD using experimental data.

IV.4.1 Experimental test case

The In this experiment, the temperature evolution of the cells is studied under various airflow conditions. Here we focus only on the results provided by the natural convection case, with no inlet velocity imposed at the inlet of the enclosure.

IV.4.2 2D simulation

IV.4.2.1 Equations and computational domain

A two-dimensional simulation of the transient conjugate heat transfer is performed to simulate the experimental case. The full set-up is represented in two dimensions, considering the horizontal plane of the real geometry. The computational domain Ω is divided into two

Parameters Values Units

Battery properties

ρ s 2227.7 kg.m -3 C p,s 837.4 J.kg -1 .K -1 λ s 0.488 W.(mK) -1

Table IV.4: Physical properties of the batteries

The source term in the solid domain Ω s represents the heat generated by the cells. Unlike the previous application case, only the term accounting for the Joule effect is considered here.

When using this type of model, it is a common assumption to neglect the term related to the entropic heating, assumed to be negligible in front of the ohmic heating. The expression of the source term in the heat equation is then:

q g ∝ R(T ) • I 2 (IV.24)
meshed solid domain As the vertical dimension is not considered in this 2D computation, the effects of gravity are not considered and the heat transfer by natural convection turns into a pure conduction problem between the air and the batteries here.

Γ 1 imposed temperature and velocity Γ 2 symmetry and adiabatic for thermal equations Γ 3 free outlet Γ 4 symmetry Γ 5 wall condition γ i coupling conditions (fluid/solid interface)

Table IV.5: Boundary conditions

IV.4.2.2 Mesh convergence for 2D simulations

First, we perform a mesh convergence study for the unstructured grids representing the fluid and solid domains. The mesh is refined using an attractor operator near the cells wall, i.e.

the solid/fluid interface, and in proximity to the plexiglass wall. The results show that TrioCFD reproduces well the predictions of Fluent. Significantly, the temperature for cells 1 and 2 (two first cells in flow direction) overlap well between the two codes. A slight difference between is to note for cells 3 and 4. The temperature difference between the two codes is slightly higher for cell 4 than for cell 3. Cell 4 is the last in the flow direction, and one can guess that there might be some difference in the velocity field prediction, leading to this temperature discrepancy.

IV.4.2.4 Deterministic temperature prediction

Here the results of the simulation in the case of pure heating and natural convection are Discussion on the source term expression A specific polynomial expression was chosen for the resistance model R(T ) in the source term q g to reproduce the experimental results.

The present reference [START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF] used a source term constructed from the fitting with their experimental measurements. This approach considers the dependency of the source term directly with the time of the experiment. First, this expression for the source was not explicitly given in the reference. Also, considering the time as input to traduce the temperature values would be equivalent to doing a backward problem and does not match the present work's objectives.

Indeed, this method is not consistent if we want to find a model of heat source representing the processes governing the actual heating of Li-ion cells. Especially, the approach described in the first chapter (see section II.2.2.2) is based on observations on the Li-ion cells behavior.

As performed in the last application and shown in literature, the resistance is often taken depending on the temperature. To use the present model from a pure design perspective, it has to consider the temperature of the cells as input, which results from the heat transfer in the immersed configuration. Now that the physical assumptions have been stated, the issue is to choose an expression of this R(T ) polynomial. Following literature considerations (as also seen in the last section), we consider a third-degree polynomial for the internal resistance:

R(T ) = α 3 T 3 + α 2 T 2 + α 1 T + α 0 (IV.25)
The choice of the source term is then reduced to a choice on the four coefficients of the polynomial expression. The construction of the polynomial expression is covered with further details in chapter VII. To assess the accuracy of the CFD model, an L2 error was computed between the transient prediction of temperature and the experimental data at some given times. Several CFD computations are performed, and the errors are compared. Then the polynomial used in the simulation presenting the smallest error was selected to conduct the simulations in this section. We denote this polynomial R CL2 (T ).

Results

First, let us focus on the resulting temperature fields during the whole simulation. We focus the attention briefly on the solid domain with a more narrow color scale, in Fig.

IV. [START_REF] Chalise | Conjugate Heat Transfer Analysis of Thermal Management of a Li-Ion Battery Pack[END_REF]. One can see that the cells 2 and 3 at the middle present higher temperatures. Because of the cell arrangement and the wall on the left side, the air around those cells is in a confined zone, and the thermal diffusion seems less important in that area. Physically, the heat exchange between the fluid and solid is less important as it is performed with temperatures close to each other around these cells. Then some hot spots appear in these cells. This fact is even more visible in cell 1, which is in contact with air from the inlet. The bottom part of the cell seems better cooled than the part facing the confined area. The temperature field inside this cell presents then a more important gradient along the vertical direction. cess. Indeed, the final temperature predicted by the model is pretty close to the experimental data. Also, the temperature given by the CFD overlaps with the measurements error envelope systematically. However, some slopes changes are occurring in the experimental data after t = 300 sec, again after t = 600 sec and at the end after t = 1400 sec. Even if the CFD stands in the error envelopes, the current model does not seem capable to capture these slope changes. The CFD solution looks like a pure heating problem with poor cooling capability.

Indeed, no airflow was generated under the current assumptions and this 2D geometry. The conjugate heat transfer turns into a pure conduction problem between the two media. In the next section, a simulation in 3D with the according hypothesis is performed to observe if the thermal behavior will change from the 2D results.

IV.4.3.1 Mesh convergence for 3D simulation

Before assessing the performances of the code in 3D, we perform a mesh convergence study for this geometry. Like in the previous section, the 3D geometry is meshed with an unstructured grid of tetrahedra. Three meshes contain respectively 6.8E + 4, 1.03E + 5 and 2.21E + 5

tetrahedra. The convergence is assessed with the temperature evolution in the probe position as in the 2D case. In the vertical direction, the probes are located at the mid-height of the batteries.

We can see the convergence is reached for those meshes. The following simulations will be performed with the mesh containing 1.03E + 5 tetrahedra. This mesh is visible through the few boundaries represented in Fig. IV.22. The simulations, which results are presented in the next section, were run on 180 cores during 6 hours, using this mesh. From what we've seen with the slight cooling process occurring with natural convection, and the resulting behavior in the far airfield, we can expect that the 2D simulation overestimates the temperature prediction compared to the 3D results. The Rayleigh number is computed using the same Prandtl number as the previous test case.

We have Ra = 1.4 • 10 4 . Given the Ra number, the flow around the cylinder is considered laminar. Then, the resulting Nusselt number is computed with

N u = 0.59 • Ra 1/4
And in this case, we have N u = 6.5. So for the first cell, the convective heat transfer is dominant compared to the conduction.

Let's perform the same computation for the cells inside the pack, which are supposed to be less cooled due to the slower air velocities in this zone. For cell 2, the temperature difference from bottom to top is estimated approximately to T = 0.1 [K] (if we consider a surface facing the inner pack). The resulting Nusselt number is lower than in the previous case but still above 1. We find N u = 4. The convective heat transfer is still dominant but with less importance.

IV.4.4 Comparison 2D/3D and experimental data

Now we compare the temperature prediction at the probe position on cell 1 between the The first observation is about the shape of the temperature predicted by the 3D CFD simulation. For the first 300 seconds, the temperature of the 2D and 3D solutions overlaps pretty well. However, as time increases, the slope of the 3D solution decreases slightly. This behavior could be explained by the convective heat transfer between the solid and fluid due to the flow generated around cell 1. The comparison with the curve of the 2D simulation, which is a straight line, supports this explanation. More precisely, we notice also that the curve of the 3D CFD presents slope changes similar to the experimental data. The slope break after t ≈ 600 sec is visible on those two curves. It seems the 3D computation allows reproducing a richer physics occurring in the experiment and taking into account the cooling effects due to natural convection.

However, on the other hand, the values of temperature predicted in 3D go out of the measurement errors envelope after t = 500 seconds. The final temperature prediction in 2D appears closer to the experimental data than in 3D. However, the previous section's plots of temperature velocity vector fields showed that the 2D simulation would overestimate the temperature prediction by omitting the natural cooling. The source term is expected to be the most responsible for the temperature evolution slope. The behavior of the model selected for the source term in the solid domain q g , through the polynomial R CL2 (T ) seems not accurate enough with respect to the experimental data if used with the 3D computation. Indeed, the R CL2 (T ) model was calibrated using the 2D configuration. Calibration with the 3D computations would be necessary to fit the temperature evolution on the experimental data.

However, the runs with the 3D geometry are quite expensive in terms of computational cost, and a proper calibration process turns out to be prohibitive with this configuration. Further discussion on the calibration of the source term, including uncertainties, and the choice of the computation geometry is given in the chapter VII.

IV.5 Chapter conclusion

In this chapter, we described the CFD approach for the simulation of Li-ion batteries immersed in air. The set of equations solved is detailed by considering the physical and modeling assumptions. The numerical methods employed in the code TrioCFD are described to give an overview of the code capabilities. The second test case is the validation of the code with experimental measurements of temperature in natural convection, using a test case from literature. The comparison between 2D and 3D simulations highlights some of the important physical processes to consider when predicting the batteries' temperature evolution. Besides, this configuration is also used to perform another comparison with Fluent on a forced convection regime. More importantly, the transient temperature predictions are compared with 2D and 3D models and experimental data. Finally, using this more straightforward physical configuration (natural convection), we assess the ability of the code to reproduce the heating of Li-ion batteries. We learn that both the effects of the pure heat transfer between the media and the way we represent the heating source of the batteries will affect the temperature prediction. Regarding the absolute value of the temperature, for instance, the maximal temperature reached at the end of the simulation, the heat source term from the solids seems to be the most important parameter.

In the following chapters, further analysis methods, including uncertainties and statistical observations, allow addressing these assumptions.

V.1 Uncertainty Quantification for computer models V.1.1 Issues and goals of Uncertainty Quantification (UQ)

So far, this manuscript has covered the main physical issues related to the cooling of Li-ion batteries and numerical approaches to address the deterministic simulation of such problems.

The objective of numerical simulation is to show that the constructed models can reproduce experimental measurements of the quantities of interest. Indeed, the experimental mea- In many engineering and real-world applications, it turns out that some of the involved parameters are pretty difficult to estimate. Often, there is a lack of experimental data to evaluate the relevant parameters. When available, the experimental measurements come with errors too. Also, sometimes the choice of the modeling approach forces us to make assumptions to represent the physical process. So, it is crucial to select values and models relevant to the physical phenomenon to be reproduced.

These considerations highlight the need of taking into account uncertainties in the simulations. As it is impossible to alleviate completely the unknowns related to some problem, the answer to these issues is to deal with the uncertainties. The challenge of Uncertainty Quantification (UQ) is then to answer the following questions: How can we incorporate the uncertainties in the physical input parameters to increase the predictive character of the simulations? What are the impacts of the uncertainties of input parameters on the predicted quantities of interest? Can we use the experimental data to learn the model's behaviour and reduce the uncertainty on the inputs?

The uncertainties may arise from many different sources and then be of different natures: model or experimental measurement errors, modeling errors due to simplifications/assumptions, inherently variable parameters [START_REF] Der | Aleatory or epistemic? Does it matter?[END_REF]. However, it is common to draw two distinct categories of uncertainties: aleatoric or epistemic uncertainty.

Aleatoric uncertainty represents the intrinsic randomness of a phenomenon. These uncertainties cannot be reduced by the knowledge brought by experts or any supplementary investigations. The studied process is random, and its outcome can only be assessed with a certain probability.

Epistemic uncertainty refers to a lack of knowledge in the investigated phenomenon. This uncertainty originates from a complex physical process not fully understood yet. The amount of unknown information can be reduced by further experimental or numerical investigations, improved conceptual models or by the experience of some experts on specific matters.

As further explained in [START_REF] Der | Aleatory or epistemic? Does it matter?[END_REF], the type of uncertainty considered depends mainly on the perspective chosen by the modelers. However, it is still relevant to define it correctly with respect to the problem to evaluate which parameters uncertainty can be reduced.

Naturally, the UQ methods rely on a mathematical framework based on probability theory and statistics to tackle these issues.

"Uncertainty Quantification is the coming together of probability and statistics with the real world.", T. J. Sullivan in [START_REF] Sullivan | Introduction to Uncertainty Quantification[END_REF].

Probability theory is the way to deal with the uncertainty related to some parameters. In the UQ framework, no parameter or simulation outcome is determined 'surely', but with an associated probability as the values of parameters vary in some range, infinite or finite. The parameters are random variables, taking their values according to a given probability density function (PDF). The quantities of interest are estimated by evaluating statistical moments such as mean, variance, quantiles and so on, which are quantifying the variability of the studied model.

V.1.2 Some objectives in Uncertainty Quantification

From a general point of view, let us assume we study a physical process represented by a

numerical model M. The model F maps a set of m inputs X = (X 1 , ..., X m ) into a set of p quantities of interest Y = (Y 1 , ..., Y p ) such that Y = M(X) (V.1)
The two main objectives addressed in this thesis are described in the following. Note that a consistent mathematical notation is provided in the next section.

Forward propagation

Let us assume that some uncertainties are associated with the vector of inputs X. Then, how much will the uncertainty on the inputs impact the outputs' variability in Y? Practically, the objective is to compute statistics of the quantities of interest in Y such as mean, variance or probability density functions, to understand the behavior of the model M under the considered uncertainties on the inputs X.

Model calibration Knowing some outputs Y true taken for granted as the 'true' output values of the physical process (for instance experimental data), and assuming the model M is able to reproduce the real data with an accepted discrepancy ε, i.e.

Y true = Y + ε (V.2)
We are interested to find out which values of inputs denoted X calib would lead to a response from the model M the closest to this data in some norms.

This chapter aims to explain the methods and tools used to perform UQ studies. For comprehensiveness, the required probabilistic tools mentioned earlier are described formally. Then each of the following sections corresponds to deeper explanations on the methods employed to tackle the UQ problems: sampling, Monte Carlo methods for forward propagation, surrogate modeling approach, global sensitivity analysis and finally Bayesian calibration.

V.2 Some mathematical definitions and UQ workflow

For further details on the definitions and properties of the mathematical objects defined in the following, the reader can refer to [START_REF] Sullivan | Introduction to Uncertainty Quantification[END_REF].

Let Ω be a set called sample space, F a σ-algebra on Ω, i.e. a collection of subsets of Ω including the empty set ∅. The elements of Ω are called measurable sets or events and the space defined by (Ω, F) is a measurable space. We define a signed measure µ as a function on the measurable space (Ω, F) such that:

µ : Ω → R ∪ {+∞ OR -∞} (V.3)
and µ(∅) = 0 µ is called a probability measure if µ(Ω) = 1 and we denote it P. Then we define the space (Ω, F, P) as a probability space.

V.2.1 Random variables and vectors

Let (Ω, F) and (X , G) measurable spaces. A function f : Ω → X generates a σ-algebra on Ω by:

σ(f ) := σ {f ∈ E|E ∈ G} f is called a measurable function if σ(f ) ⊆ F.
If we define such a function X on a probability space (Ω, F, P), and if we take X ⊂ R, X is a real valued random variable defined by:

X : Ω → X ω → X(ω) (V.4)
We note x = X(ω) a realization of the random variable X from an event ω ∈ Ω. X induces a probability space (X , B(X ), P X ) where B(X ) is a Borel σ-algebra of the event space X . The probability measure P X on this space is defined by:

P X (ω ∈ E) := P(X ∈ E), for E ∈ G (V.5)
P X is the distribution or law of the random variable X. It is defined by a so-called Cumulative Density Function (CDF) F X or Probability Density Function (PDF) π X given by the following expression, for a realization x:

F X (x) = P X (X ≤ x) (V.6) π X (x) = dF X (x) dx (V.7)
In the following, we suppose that the output of random variables have a finite variance and belong to the Hilbert space of square integrable functions.

The PDF π X is used to get the probability of some events I ∈ X by computing:

P X (X ∈ I) = I π X (x)dx (V.8)
It can be seen as the continuous sum of the realizations x on the subset I defining an event of X . Similarly, the statistical moments associated to the random variable X are defined as follows. The first order moment is the mean. In the following we denote it by the operator E and it is given by:

E[X] = X x • π X (x)dx (V.9)
The second order moment is the variance written with the operator V and computed by:

V[X] = X x 2 • π X (x)dx (V.10)
In the field of statistics it is common to use the information provided by the variance through the so-called standard deviation defined as:

σ X = V[X] (V.11)
It is useful for the analysis of any process as its physical dimension is the same as the mean and the outputs of the random variables itself.

Random vectors

In a similar fashion, we can define random vectors. They are vectors containing a finite collection of random variables.

X : Ω → X ⊆ R m (V.12) ω → X(ω) = (X 1 (ω), ..., X m (ω))
where each X i is a random variable. The Cumulative Density Function for the random vector can be defined by:

F X = P(X 1 ≤ x 1 , ..., X m ≤ x m ) (V.13)
where the x i are realizations of each random variables X i .

The first and second order statistics coming from the random vectors are given by the following:

E[X] = X x • π X (x)dx (V.14)
where x = (x 1 , ..., x m ) is the vector containing the realizations of the random variables

X 1 , ..., X m .
The second order moment is the variance written with the operator V and computed by:

V[X] = X x 2 • π X (x)dx (V.15)
The joint probability density function π X of the random vector is defined by:

π X = π x 1 ,...,xm (x 1 , ..., x m ) (V.16)
In the following of the manuscript, the bold capital letters will designate random vectors (X).

The random variables will be designated by normal font capital letters (X). The realization of random vectors is designed by small bold letters (x) and normal font for random variables (x).

Covariance For a set of two random variables X 1 and X 2 , the covariance measures the simultaneous variation of the two random variables. The covariance is defined by:

Cov(X 1 , X 2 ) = E[(X 1 -E[X 1 ])(X 2 -E[X 2 ])] (V.17)
Note that, by definition we have Cov

(X 1 , X 1 ) = V[X 1 ]
. This notion is generalized to a larger set of m random variables X = (X 1 , ..., X m ), using the covariance matrix K. The covariance matrix is used to represent the variance of the full vector of random variables. The expression of the matrix is:

K =         V (X 1 ) Cov (X 1 , X 2 ) • • • Cov (X 1 , X m ) Cov (X 2 , X 1 ) . . . • • • . . . . . . . . . . . . . . . Cov (X m , X 1 ) • • • • • • V (X m )         (V.18)
It is a symmetric matrix so K T = K and the diagonal contains the variance of each random variable.

Correlation matrix

The correlation between two variables X 1 and X 2 is given by:

Cor(X 1 , X 2 ) = Cov(X 1 , X 2 ) σ 1 σ 2 (V.19)
We define the correlation matrix for m random variables X = (X 1 , ..., X m ) as:

R =         1 Cor (X 1 , X 2 ) • • • Cor (X 1 , X m ) Cor (X 2 , X 1 ) . . . • • • . . . . . . . . . . . . . . . Cor (X m , X 1 ) • • • • • • 1         (V.20)

V.2.2 Workflow of a typical UQ problem

Considering the previous defined mathematical tools, let us describe more in detail the standard steps to go through when solving a UQ problem [START_REF] Zhang | Basic Framework and Main Methods of Uncertainty Quantification[END_REF]. The workflow is summarized in Fig. V.2 with the relevant notations described in this chapter. Note that what follows is only one of many different ways of sketching the different UQ steps in the literature. From a numerical perspective, the software UQLab has been used in the framework of this thesis, to perform all the UQ steps described here [START_REF] Uqlab | from the Chair of Risk, Safety and Uncertainty Quantification of ETH Zurich[END_REF].

Step 0: Construction of the numerical model The physical process is usually approximated with a numerical model M, solving a set of discretized equations. This model is seen as representing reality with a discrepancy from the real process. This discrepancy can be explained by experimental measurements uncertainties, model errors, imprecise calibration, numerical errors ...

Step 1: Uncertainty characterization A fundamental step is to identify and characterize system uncertainties. This step is problem-dependent and requires expertise both in the modeling and experimental aspects. More specifically, it is crucial to identify the structure of the inherent uncertainties in terms of distribution, potential correlation, and confidence intervals. Note that this analysis impacts both inputs and outputs of the numerical model.

This aspect can seem contradictory with respect to what is indicated in Fig. As explained earlier, uncertain inputs are then represented by random variables. Uncertainties will be characterized by some probability density functions, relying on specialists in the related field. The choice of the so-called prior distribution is expected to impact the variability of the output.

Step 2: Forward propagation of the uncertainties Once the uncertainties on the inputs are formally represented through the random variables and their associated probability density functions, the problem is to compute statistics on the output variability. Given an uncertain input, represented by the random vector X = (X 1 , ..., X m ), we consider the scalar quantity of interest (QOI) given by a random variable Y ∈ Y (resp. vector of QOI) such that Y = M(X). The forward propagation consists in computing the statistics, for instance, mean and variance, of the output Y .

There are two main families of approaches denoted as intrusive and non-intrusive methods. The intrusive methods consist in inserting random variables inside the model M and making its output random [START_REF] Kusch | Intrusive methods in uncertainty quantification and their connection to kinetic theory[END_REF]. This methodology requires making supplementary implementations in the numerical model, as stochastic processes represent the unknown of the system equations. Non-intrusive techniques consider the model M as a black box. Specifically, the computation of the model is seen as deterministic: for a given input, the output will always be the same. The variability comes from the uncertainty of the input values. The nonintrusive techniques present a lot of advantages in UQ problems with numerical models [START_REF] Francesco | Predictive numerical simulations for rebuilding freestream conditions in atmospheric entry flows[END_REF].

The model can be used as it is, regardless of its internal complexity, to perform a UQ analysis.

However, the non-intrusive approach can increase the computational resources considerably because of the so-called curse of dimensionality [START_REF] David | High-Dimensional Data Analysis: The Curses and Blessings of Dimensionality[END_REF]. In this work, only non-intrusive uncertainty quantification techniques have been performed. In V.3, we detail sampling methods and especially the most common one based on Monte Carlo (MC) sampling [START_REF] Kroese | Handbook for Monte Carlo methods[END_REF].

The mathematical methods for steps 1 and 2 used in this work are described in V.3.

For practical reasons, essentially related to computational costs (as one sample is an actual evaluation of the computational model), the model M can be replaced by an accurate and fast to compute surrogate model [START_REF] Sudret | Surrogate models for uncertainty quantification: An overview[END_REF]. The construction of such model is covered in V.4.

Step 3: Sensitivity analysis This step provides further insights on the behavior of the model M. The sensitivity analysis objectives are to assess the impact of some input on the variability of the model's output.

Among many approaches possible, we focus here our attention on global sensitivity analysis [START_REF] Pianosi | Sensitivity analysis of environmental models: A systematic review with practical workflow[END_REF]. The global approach assesses the variability of the outputs on the whole inputs space. As the PDF and the above variance of the outputs are observed, it is possible to draw a hierarchy of which input X i has the most significant impact on the model's output variability [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]. This hierarchy can be established by using Sobol Sensitivity Analysis [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF]. The methodology to tackle this step is described in Section V.5.

Step 4: Inverse problem and calibration Once uncertainties on the inputs are modeled and related statistics on the output are known, a question is associated with validating the numerical solution. The numerical error coming from the forward propagation can be directly compared with experimental measurements and their error bars. Still, using experimental data directly could permit to infer some input parameters and improve the prediction capability of the solver.

In this context, a common calibration problem would be to assess which values of the input parameters reduce the error between the model and the data by minimizing a L2 error nom of the discrepancy. In the studies presented here, a Bayesian calibration approach is followed [START_REF] Ralph | Bayesian Calibration -What, Why And How[END_REF]. The mathematical and numerical tools to perform such calibration are described in Section V.6.

Uncertainties defined in

Step 1 are called the prior uncertainties, and the calibration outcome is the posterior uncertainties.

Step 5: Forward propagation of the posterior distributions As the posterior distributions are provided from the Bayesian calibration process, step 3 can be repeated using the learned distributions. This way, we can close the loop of the uncertainty quantification problem, as we can expect the output variance to be reduced using the posterior distributions. This overall process stands as an uncertainty reduction approach for the whole problem.

of output values M(X (k) ) k=1,...,N

3. Compute the PDF of the resulting output quantities. The PDF can be approximated by Kernel density estimation methods using the output samples M(X (k) ) [41].

4. Then the statistics are approximated using:

E[M(X)] ≈ 1 N N k=1 M(X (k) ) (V.23) V[M(X)] ≈ 1 N N k=1 M(X (k) ) -E[M(X)] 2 (V.24)
The question is now to determine how to choose the number of samples N for the inputs X to reach a proper convergence of these quantities. Efficiency of the sampling can be characterized by some criteria such as the space-filling and the convergence rate [START_REF] Barbillon | Méthodes d'interpolation à noyaux pour l'approximation de fonctions type boîte noire coûteuses[END_REF][START_REF] Thenon | Utilisation de méta-modèles multi-fidélité pour l'optimisation de la production des réservoirs[END_REF]. An optimal sampling should explore at best the space with a minimum amount of samples. In this work, two sampling approaches have been used following different perspectives.

Monte Carlo sampling

The most classical approach is to use the so-called Monte Carlo method. It is a robust approach because its convergence rate is independent of the input dimension m. This approach consists in sampling N M C values of the variables X i following their respective PDFs in a purely random fashion. This method presents no bias in the sampling process. However, the convergence is slow and requires a considerable number of samples, following a O(N -1/2 ) rate [START_REF] Caflisch | Monte Carlo and quasi-Monte Carlo methods[END_REF].

Latin Hypercube Sampling Optimal sampling techniques have been developed in the literature [START_REF] Dick | Digital Nets and Sequences[END_REF] to accelerate the convergence rate. Here, we focus on Latin Hypercube Sampling (LHS) [START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF], but many others possibilities exist, such as for example Sobol sequences [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF],

Halton sequences [START_REF] Halton | Algorithm 247: Radical-inverse quasi-random point sequence[END_REF], among others.

The LHS method works as follows. Consider a random vector X = (X 1 , ..., X m ) taking its values in the product space

I = m i=1 [a i , b i ].
To sample N points from these variables, the principle is to divide each of the intervals [a i , b i ] in N sub-intervals of equal probability. Note that the width of each sub-interval may vary in each i th direction depending on the PDF π X i associated to each variable X i . Then, within each sub-interval, one single sample is drawn randomly according to the PDF π X i as it would be done in MC sampling. The advantage of the method is to ensure that each component of the random vector is sampled in a fully stratified way, ensuring a good space exploration, independently of the dimension m of the inputs [START_REF] Mckay | A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code[END_REF].

123

Let us compare for a couple of random variables X = (X 1 , X 2 ) how the LHS and MC present different space-filling performances for the same number of samples. The two variables are following a uniform law in [0, 1] 2 : X ∼ U([0, 1] 2 ). We draw N = 20 samples of the random vector X. For each sample X (k) we plot the resulting point (X

(k) 1 , X (k)
2 ). The resulting points using the LHS sampled are the black crosses and their corresponding sub-intervals with dotted black lines in both directions. The points with the MC method are the red crosses.

It is visible in

V.4 Kriging surrogate model

The sampling method requires numerous evaluations of the model M to estimate the statistics of the outputs. However, in real-life applications and especially in CFD-based applications, M is usually a numerical model expensive to evaluate in terms of computational costs.

The evaluations of the statistics turn out to be prohibitive if we need to run thousands of simulations with the model. A possibility is then to replace the model M with a so-called surrogate model (or metamodel), which is a mathematical representation permitting to predict M. This thesis focuses on a specific methodology based on the so-called Kriging method.

Note that others interpolation and regression methods exist in the literature, as detailed in [START_REF] Sudret | Surrogate models for uncertainty quantification: An overview[END_REF]: Polynomial Chaos Expansion, Support Vector Machine etc... We use Kriging since one of the most used techniques to build surrogates efficiently by requiring a limited number of samples.

In what follows, we describe the necessary mathematical features to describe the surrogate modeling technique known as Kriging. First, we recall a formal definition of random processes. Then, details on Gaussian processes are stated. Finally, we describe the necessary elements to build Kriging surrogate models.

V.4.1 Random process

A random process, or stochastic process, is a collection of random variables defined in the same probability space (Ω, F, P). In particular, a continuous random process is indexed by a variable x ∈ X ⊆ R m . We then define a scalar random process Y as:

Y : X × Ω → Y ⊆ R (x, ω) → Y (x, w) (V.25)
For a ω 0 ∈ Ω, a realization of the stochastic process is every function x → Y (x, ω 0 ) associated to an event ω 0 ∈ Ω. In We define similarly a vector random process Y as the application mapping the inputs (x, ω)

into a multi-dimensional space R p such that Y = (Y 1 (x, ω), ..., Y p (x, ω)).

V.4.2 Gaussian Process

Gaussian processes (GP) are a particular type of random process. A Gaussian process is a collection of random variables {Y (x)} x∈R m indexed by x, where each variable follows a Gaussian distribution. For a finite collection of p random variables Y (x 1 ), ..., Y (x p ), we write:

(Y (x 1 ), ..., Y (x p )) ∼ N (µ, K) (V.26)
and in particular:

∀i = 1, ..., p , Y (x i ) ∼ N (µ i , σ 2 i )
Where µ = (µ 1 , ..., µ p ) is the vector containing the means of each random process. K is the covariance matrix associated to the Gaussian process such that:

K =         σ 2 1 Cov (Y (x 1 ), Y (x 2 )) • • • Cov (Y (x 1 ), Y (x p )) Cov (Y (x 2 ), Y (x 1 )) . . . • • • . . . . . . . . . . . . . . . Cov (Y (x p ), Y (x 1 )) • • • • • • σ 2 p         (V.27)
From a practical point of view, the Gaussian character of the generated random variables coming from the stochastic process can be visible in the example illustrated in Fig. V.5. We consider the one dimensional Gaussian process Y (x) indexed by a variable x ∈ [-4, 4] ⊆ R.

For a fixed index x 0 , represented by the vertical dark line, points formed by the intersection of the realizations of GP and the line x = x 0 (left) are distributed following a Gaussian distribution (right), which the mean is given by µ 0 = Y (x 0 ) and the variance by σ 2 0 . As the number of realizations increases (the number of MC samples), the sampled Gaussian distribution becomes more accurate and representative of a 'true' Gaussian distribution.

-4 -3 -2 -1 0 1 2 3 4 x -4 -3 -2 -1 0 1 2 3 4 Y (x

V.4.3 Gaussian process regression -Kriging

The use of GPs for interpolation/regression and classification problems is comprehensively introduced in [START_REF] Rasmussen | Gaussian Process for Machine Learning[END_REF]. Two properties make the Gaussian distributions and then GPs useful for regression and interpolation of mathematical models:

1. For a Gaussian random process X = (X 1 , ..., X n ), we have for a matrix A ∈ R p×n :

X ∼ N n (µ, K) ⇐⇒ AX ∼ N p (Aµ, AKA T )
It means that in practice, the sum of Gaussians is still Gaussian, and the marginal distributions of multivariate Gaussians are Gaussian.

2. The conditional distributions between two Gaussian variables are still Gaussian.

The Kriging technique was first developed by Krige in [START_REF] Krige | A statistical approach to some mine valuation and allied problems on the Witwatersrand[END_REF] as a spatial interpolation tool.

The theory behind this tool was formalized by [97]. It consists of a stochastic algorithm that relies on Gaussian processes, giving the unknown evaluation of a function at some location from already known estimations of the function [START_REF] Ginsbourger | Multiples métamodèles pour l'approximation et l'optimisation de fonctions numériques multivariables[END_REF].

In the following, we describe the key elements constituting the Kriging models, inspired from the full explanations given in [START_REF] Lataniotis | UQLab user manual -Kriging (Gaussian process modeling)[END_REF]. As the objective here is to construct a model being able to replace the true numerical model M, the Kriging model is denoted M K . Note that the Kriging models developed in this thesis considered noise-free data, and then the models act as interpolation functions rather than regression functions. The Kriging computes an outcome of the model M K , from a vector input x ∈ R m by the following expression:

M K (x) = β T f(x) + σ 2 Z(x, ω) (V.28)
The terms of the Kriging expression decompose as:

• β T f(x) is the mean of the Kriging, also called the trend.

β = (β 1 , ..., β b ) is a vector of b constant coefficients. f = (f 1 , ..., f b ) is a set of arbitrary functions.
• σ 2 is the constant variance.

• Z(x, ω) is a stationary Gaussian process, indexed by x. A stationary stochastic process means that the random variables Z(x, ω) and Z(x , ω), indexed in x and x have the same distributions.

• The probability space (Ω, F, P) is represented through the parameter ω and defined by the correlation function R defined as R : (x, x ) → R(x, x ; θ) . θ is a vector containing some of the so-called hyperparameters defined just later. This terminology is widely employed in the framework of approximating complex computer models. This set is defined as follows, for a number N DOE of observed data:

T DOE = {(X DOE , Y DOE )} ⊆ R N DOE × R N DOE (V.29)
More precisely , the DOE T DOE contains the pairs composed by the input and the corresponding value of the true model, for each observation:

X DOE = x (i) , i = 1, ..., N DOE (V.30) Y DOE = M(x (i) ) , i = 1, ..., N DOE (V.31)
From the assumptions provided by the Gaussian processes properties, the vector formed by the prediction M K (x) on the new input x and the set of observed values with the model Y DOE is Gaussian. The corresponding mean and variance are given by:

M K (x) Y DOE ∼ N N DOE +1 f T (x)β Fβ , σ 2 1 r T (x) r(x) R (V.32)
The newly introduced terms are defined as follows:

• F is a matrix containing the values of the functions defining the trend at the observed points:

F i,j = f j (x (i) ) i=1,...,N DOE ,j=1,...,b
• r(x) is the vector containing the evaluations of the correlation function between x and all the construction points:

r(x) = r i = R(x, x (i) ; θ) i=1,...,N DOE
• R is the correlation matrix defined by: R i,j = R(x (i) , x (j) ; θ) i,j=1,..

.,N DOE

The expression of the mean and variance of random variable defining the predictor M K reads:

µ M K (x) = f(x) T β + r(x) T R -1 (Y DOE -F β) (V.33) σ 2 M K (x) = σ 2 1 -r T (x)R -1 r(x) + u T (x) F T R -1 F -1 u(x) with β = F T R -1 F -1 F T R -1 Y DOE and u(x) = F T R -1 r(x) -f(x).
Finally, the Kriging predictor in x is a Gaussian random variable defined by:

M K (x) ∼ N µ M K (x), σ 2 M K (x) (V.34)
From a practical point of view, evaluating the Kriging model M K on a new point is reduced to compute the expressions in Eq. (V.34). It just consists of some matrix-vector product computations. From this perspective, one can easily see why Kriging is a convenient tool for surrogate modelling to replace the expensive model M. Now the question is to ensure that the Kriging model can reproduce well the results of the true model, based on the observations in T DOE . In other words, how are determined the hyperparameters β, θ and σ 2 involved in the above expressions and given the choice made for the functions f.

V.4.3.1 Trend of the Kriging

First, the choice to be made is about the trend of the Kriging model, i.e. the parameters in β and the form of the functions in f. Three main approaches are considered for the trend, defining in the mean time the type of Kriging used:

Simple Kriging Here the trend is assumed known with constant parameters and known functions. The parameters are all equal and β i = 1, i = 1, ..., b. The trend reads then:

β T f(x) = b i=1 f i (x) (V.35)
Ordinary Kriging Here, the trend is constant, so β i = β 0 and f i = 1 for i = 1, ..., P . But it's needed to determine the value of the parameter β 0 .

β T f(x) = β 0 (V.36)
Universal Kriging This is the most versatile and general formulation. The values of the coefficient are to determine, and the functions are usually polynomials which coefficients are also to be found.

β T f(x) = b i=1 β i f i (x) (V.37)
The trend is called linear when the functions f i are linear, e.g. first, degree polynomials for each input x i .

V.4.3.2 Correlation functions

The next really important feature of the Kriging model to define is the correlation function R : → R(x, x ; θ) leading to the construction of the correlation matrix R. The correlation function will determine how much the model should consider the close presence of known data points in predicting the new outcome in x. In other words, it quantifies the similarity between observations and new points, depending on the distance between them. They are expressed as functions of the distance between the two considered points: |xx |. Also, to get the best model possible, the behavior of the model for two close points x, x needs to be traduced by the behavior of the covariance function [START_REF] Thenon | Utilisation de méta-modèles multi-fidélité pour l'optimisation de la production des réservoirs[END_REF]. Many covariance functions are used in literature. We write here the expression for the most common ones: linear, Matérn, exponential and Gaussian.

Note that the correlation functions need to respect some conditions to be suitable for Kriging: the function is symmetric, and the correlation matrix induced by the evaluations of the function is positive semi-definite. For the sake of simplicity and illustration, we consider in this specific section a model taking a scalar x as input: x is no longer a vector in R m , and we write x = x ∈ R. Also, the vector of hyperparameters θ becomes a scalar θ in this specific case. Note that we can define a correlation function for vector input variables, as detailed in the last paragraph of this section.

So we consider two distinct points x, x ∈ R and a hyperparameter θ. The expressions of the correlation functions are given in Tab. V.1:

Linear R(x, x ; θ) = max 0, 1 -|x-x | θ Exponential R (x, x ; θ) = exp -|x-x | θ Gaussian R (x, x ; θ) = exp -1 2 |x-x | θ 2 Matérn 3/2 R(x, x ; θ) = 1 + √ 3 |x-x | θ exp - √ 3 |x-x | θ Matérn 5/2 R (x, x ; θ) = 1 + √ 5 |x-x | θ + 5 3 |x-x | θ 2 exp - √ 5 |x-x | θ Table V.1: Expression of correlation functions.
The plot of these functions against the distance |x -x | can be found in many references such as [START_REF] Rasmussen | Gaussian Process for Machine Learning[END_REF][START_REF] Lataniotis | UQLab user manual -Kriging (Gaussian process modeling)[END_REF]. The hyperparameter θ is often called the length scale. Relying on the expressions of the above functions, these parameters can be seen as the quantification of how much the response of the process should vary from one point to the other. If θ is small, the correlation functions present strong gradients in the |x -x | direction, and the Gaussian process might vary a lot between those two indexes. On the opposite, a large θ will produce a smooth Gaussian process. An example from [START_REF] Lataniotis | UQLab user manual -Kriging (Gaussian process modeling)[END_REF] is shown in Fig. V.6 with a Gaussian correlation function.

For processes with vector inputs, we can similarly define multi-dimensional correlation functions. Then the vector of hyperparameters θ contains the lengthscales θ i for each function.

Two families of function are considered: separable or ellipsoidal correlation functions. This estimation method is based on the maximization of the likelihood function. This function is dependent on the parameters of a statistical model (here the hyperparameters) and computed from some observed data (here the data in T DOE ). In other words, for a random variable X following a PDF parametrized by a scalar u, the likelihood L is a function of u given an observation of the process x = X, and it is written:

L(u|x)
It gives the probability of observing x if the parameter was u.

In the case of the Kriging, the expression of the likelihood is given in Eq. (V.40). It is expressed in function of the hyperparameters: β, σ 2 , θ and the observed data Y DOE .

L β, σ 2 , θ;

Y DOE = (det R) -1/2 (2πσ 2 ) N/2 exp - 1 2σ 2 (Y DOE -Fβ) T R -1 (Y DOE -Fβ) (V.40)
The expression of the vector parameters maximizing the likelihood function, β and σ are written in Eq. (V.41). Full details on the calculation to obtain these expressions are given in [START_REF] Thomas | The Design and Analysis of Computer Experiments[END_REF]:

β = β(θ) = F T R -1 F -1 F T R -1 Y DOE σ 2 = σ 2 (θ) = 1 N (Y DOE -F β) T R -1 (Y DOE -F β). (V.41)
These expressions still show the dependency on the length scales of the correlation function θ. The hyperparameters θ are obtained by minimizing the log of the Likelihood:

θ = arg min θ∈D θ [-log L(θ; Y)] (V.42)
Constructing the Kriging model is then reduced to solve this optimization problem.

V.4.3.4 Quality assessment

Here, we describe the common tools used to assess the quality of the constructed Kriging. In other words, how well it is capable of fitting the true model M.

Validation error

This technique is used when other data are available aside from the DOE.

This other set of data is called validation set denoted T val . The validation set is constructed as the DOE: we sample N val inputs X val = x (i)

val , i = 1, ..., N val . Then we evaluate the true model M on these inputs. The supplementary computational cost is N val evaluations of the true model. Once the Kriging model M K is constructed, the idea is to evaluate it on the inputs from the validation set X val . Several quantities can be computed from this analysis with the validation plot. For instance the error ε val between the values computed by the Kriging and the true models can be estimated using the coefficient of determination [START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF] such as: [START_REF] Vahid Esfahanian | Capacity Fade Analysis of Lithium-Ion Cells Utilizing Uncertainty Quantification Approach[END_REF] where m is the mean of the true model responses on the validation set: m = 1

ε val = 1 - N val i=1 M(x (i) val ) -M K (x (i) val ) 2 N val i=1 m -M(x (i) val ) 2 (V.
N val N val i=1 M(x (i) val ).
A value close to one means that the Kriging model can reproduce the data from the true model in the validation set.

The validation set has to be chosen with caution to avoid bias in the validation procedure.

A good practice is to build T val using an independent sampling set. As in practice, the evaluations of the true model are expensive, a small sample set using a space-filling technique such as LHS is constructed. However, the risk is that some sampled points could be too close to the ones used in the construction, and the overall error might be underestimated. A way around to enforce the validation is also to compute a cross-validation error, and more particularly the one called Leave-one-out error [START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF].

Leave-one-out error

This technique is a handful as it does not require new evaluations of the true model. Based on the initial DOE, the idea is to compute the error between the true model values and a Kriging model, which was constructed removing one point in the DOE.

The Leave-one-out error ε LOO is computed as:

ε LOO = 1 N DOE • N DOE i=1 M(x (i) ) -M K {-i} (x (i) ) 2 V[Y DOE ] (V.44)
where V[Y DOE ] is the statistical variance of the model observations and M K {-i} the Kriging model constructed on a DOE deprived of the point x (i) . The drawback of this technique is that it is required to construct as many Kriging models as points in the DOE, which could be costly in some cases (see section V.4.3.5).

Physically, for each input x (i) , we sum up the deviations between response from the true model and a Kriging model constructed without this input. If the Kriging is well constructed, it is expected that M K {-i} (x (i) ) is not too far from M(x (i) ) and the error computation is a sum of terms close to zero. Then we want the overall error ε LOO to be as close to zero as possible.

V.4.3.5 A word on multiple outputs Kriging

So far, we have discussed the theory of Kriging models with output scalar values. However, in the physical applications covered in this thesis, the quantities of interest are time-dependent.

As seen in the chapters showing the performance of the numerical models, the quantities outputted by the code are the temperature evolution. A common practice to represent time dependant processes with Kriging models is to consider the output of the model M K as a vector containing the evaluations of the quantity of interest at different times. We consider here a numerical model which output vector values such that:

M : R m → R p X → Y = M(X) (V.45)
The overall process for the Kriging construction stays the same. However, the predictors given by the Kriging model are no more scalar values but vectors. Then, each predictor is defined as for a single output Kriging, independently of the others.

It is required to have a set of p vectors for each type of hyperparameter.

µ M K (x) = (µ 1 (x), ..., µ p (x)) (V.46) σ 2 M K (x) = σ 2 1 (x), ..., σ 2 p (x) (V.47)
where the µ i (x) and σ 2 i (x) are the mean and variance of each predictor as defined in Eq. (V.34). Therefore, the construction of the Kriging with p outputs is made by solving p optimization problems as in Eq. (V.42) to get values of each hyperparameter vector.

V.5 Global sensitivity analysis

The global sensitivity analysis aims to quantify the effect of each uncertain input and their correlated effect on a given scalar quantity of interest computed by the model [START_REF] Chen | Analytical Variance-Based Global Sensitivity Analysis in Simulation-Based Design Under Uncertainty[END_REF]. The analysis is performed using a variance-based decomposition of the black box model known as Analysis of Variance (ANOVA) decomposition [START_REF] Sullivan | Introduction to Uncertainty Quantification[END_REF]. Here, we consider the model M giving a single scalar y as output, from the evaluation of the m inputs x = (x 1 , ..., x m ). The ANOVA decomposition of the model M evaluated in a realization x can be written as follows:

M(x) ANOVA := M 0 mean + m i=1 M i (x i ) f irst order + m i 1 =1 m i 2 =i 1 +1 M i 1 i 2 (x i 1 , x i 2 ) second order + • • • + M 1,...,m (x 1 , ..., x m ) m-th order . (V.48)
The numerical model is seen as a sum of functions, representing the contribution from each input variable x i (the functions M i ) and their correlations (M ij etc) up to the m-th order. This decomposition allows to perform a sensitivity analysis based on the evaluation of the so-called Sobol indices [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF]. Each input is seen as a realization of the random vector X = (X 1 , ..., X m ). The Sobol indices compute the ratio of the conditional variance over the variance of full model response. More specifically, the first order Sobol index S 1 i for an input variable X i is defined by:

S 1 i = V [E[M(X)|X i ]] V [E[M(X)]] (V.49)
The first order index quantifies the effect of the variable X i only on the variability of the quantity of interest Y by computing the variance ratio in Eq. (V.49). The conditional variance E[M(X)|X i ] express the variability of the mean of the model response conditioned to the parameter X i .

The total Sobol index S T i for the variable X i is computed as follows:

S T i = 1 - V [E[M(X)|X ∼i ]] V [E[M(X)]] (V.50)
where X ∼i represents the set of all input variables and their combinations except X i . The total Sobol quantifies the impact of the variable X i and its correlation with all other variables on the variability of Y .

For a given input X i , a Sobol index close to 1 means that this variable significantly impacts the variability of the quantity of interest Y . On the opposite, a Sobol index close to 0 means that the quantity Y is almost independent of the variability of X i . Note that for a given quantity of interest, all the Sobol indices sum to 1.

To compute those indices, here are the expression of the involved variances. The total variance of the model output V [E[M(X)]] is computed using:

V [E[M(X)]] = X M(x)π(x)dx -M 0 2 (V.51)
This expression involves the multidimensional integration of the product between the model

M and PDFs π(x).

The conditional variance for the first order index are given by:

V [E[M(X)|X i ]] = X ∼i M(x)π(x ∼i )dx ∼i -M 0 2 (V.52)
where the notation ∼ i designates all the input variables except the i th one:

x ∼i = (x 1 , ..., x i-1 , x i+1 , ..., x m )
Note that for vector outputs of the model, i.e. the vector quantities of interest Y = M(X) ∈ R p , total and partial Sobol index are defined for each of the scalar quantities of interest Y j composing Y.

In practical applications, the MC methods seen in the forward propagation section V. Note that the uncertain inputs X of the model and the discrepancy ε are treated independently. So we have for the PDFs of both random vectors:

π(x I ) = π(x, x ε ) = π(x)π(x ε )
Then the expression of the posterior distribution is given by:

π(x, x ε |y exp ) = 1 Z π(x)π(x ε )L(x, x ε ; y exp ) (V.56)
This expression depends on the likelihood function which is computed by:

L(x, x ε ; y exp ) = 1 (2π) p det(K -1 (x ε )) exp - 1 2 (y exp -M(x)) T K -1 (x ε ) (y exp -M(x)) (V.57)

V.6.2 Solution and posterior predictions

One can compute the so-called point estimates to summarize the posterior distribution. The statistics usually observed are the mean of the posterior distributions and the maximum a posteriori (MAP). The mean of the distribution is expressed by:

E(X I |Y exp ) = X I x I π(x I |y exp )dx I (V.58)
The MAP gives the mode of the posterior distribution. In short, it gives the values of inputs parameters, denoted x M AP I , maximizing the likelihood function:

x M AP I = arg max x I L(x I ; y exp ) (V.59)
In practical problems, the posterior distribution is used to compute the expectation of the quantities of interest given by the model. Let's consider a scalar QOI y j ∈ y = M(x) for j = 1, ..., p, and note M j the evaluation of the model corresponding to this sole QOI. Then we can evaluate the expectation of the model response from the posterior distribution:

E[M j |y exp ] = X I M j (x I )π(x I |y exp )dx I (V.60)
Also, the amount of uncertainty related to the posterior estimation of the input parameters can be estimated with the covariance matrix of the random vector X, or with the variance of each scalar input parameter X j . The variance of such random variable is expressed by:

V(X j |Y exp ) = X j x 2 j π(x j |y exp )dx j (V.61)
By computing the ratio of variances of the prior distribution π(X I ) and the posterior, one can assess the amount of uncertainty reduction that has been made with the Bayesian calibration.

Such results are given for the application in chapter VII.

Finally, the process of forward propagation described in V.3 can be repeated by integrating along the posterior distribution. Then the statistics such as mean and variance of the model predictions are given with a sampling made following the posterior distribution. The distributions of the model M outputs can be computed using both prior and posterior distributions.

V.6.3 Numerical resolution: MCMC algorithms

In real applications problems, the computation of the posterior distribution as described in Eq. (V.56) is intractable in practice. The way around is to use sampling methods. Based on the Glivenko-Cantelli theorem, one can approximate the (unknown analytically) posterior π(x I |y exp ) by sampling a large numbers T of realizations x t I for t = 1, ..., T [START_REF] Donnet | Bayesian Inference for Latent Variable Models. en. Online course[END_REF]. The distribution of these sampled realizations converges towards the desired posterior distribution.

V.6.3.1 Markov Chain Monte Carlo

The most common approach is to use Markov Chain Monte Carlo algorithms to simulate this desired posterior distribution. This method constructs a Markov chain of the inputs realizations (x 1 I , x 2 I , ...), with an invariant distribution being the sought posterior distribution. Markov chains are defined by their transition probability to move from the step t to t + 1, denoted:

K(x (t+1) I |x (t) I )
The invariant distribution of the chain is the posterior one if the Markov Chain is reversible, i.e. if it fulfills the condition:

π(x (t) I |y exp )K(x (t+1) I |x (t) I ) = π(x (t+1) I |y exp )K(x (t) I |x (t+1) I ) (V.62)
The probability to move from the iteration t to t + 1 is equal to the probability to move back from t + 1 to t. Then, the distribution at iteration t + 1 is given by:

π(x (t+1) I |y exp ) = X I π(x (t) I |y exp )K(x (t+1) I |x (t) I )dx (t) I (V.63)

V.6.3.2 Metropolis-Hastings and Adaptive Metropolis algorithms

Several algorithms allow to construct Markov Chains with this property: Metropolis-Hastings,

Adaptive Metropolis [START_REF] Heikki | An Adaptive Metropolis Algorithm[END_REF], Hamiltonian Monte Carlo [START_REF] Nagel | Hamiltonian Monte Carlo and Borrowing Strength in Hierarchical Inverse Problems[END_REF], Affine invariant ensemble algorithm [START_REF] Goodman | Ensemble samplers with affine invariance[END_REF], among others. A comprehensive list is available at [START_REF]MCMC algorithms[END_REF]. The most classical is the Metropolis-Hastings algorithm (MH) [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov Chains and their applications[END_REF]. The other algorithms rely on the principle of the MH algorithm with some enhancements. In the applications of this thesis, the MCMC algorithm used was the Adaptive Metropolis. Here we describe first the MH algorithm and then detail the specificities brought by the AM approach.

Metropolis-Hastings

The basic principle of the method is to explore the posterior distribution with a random walk of T iterations. Each iteration consists in drawing a candidate x (c) I using a proposal distribution to move from iteration t to t + 1 [START_REF] Donnet | Bayesian Inference for Latent Variable Models. en. Online course[END_REF]. The new candidate is accepted with a probability computed based on the proposal distribution. The operation of the algorithm is given in Alg. Now, we can approximate the model response for the QOI y j from the posterior distribution.

That is the integral in Eq. (V.60), by using the T samples of the Markov Chain:

E[M j |y exp ] ≈ 1 T T k=1 M j (x (k) I ) (V.65)
The difficulty in this algorithm is to choose properly the proposal distribution π(x

(c) I |x (t)
I ). The issues related to the choice of the proposal distribution is discussed in the next paragraph, allowing us to introduce the advantages of the Adaptive Metropolis algorithm.

Adaptive Metropolis

A good proposal distribution is close to the searched posterior distribution. However, in practical problems, we do not have any expectations of the shape of this distribution.

A common approach is to choose a Gaussian distribution for the proposal, centered and with a covariance matrix K. If the variance is too large, the next candidate will often be rejected, and the chain will not move. Conversely, a too-small variance could lead to too many accepted candidates. The chain can be stuck in a region and not explore well the support of the posterior [START_REF] Donnet | Bayesian Inference for Latent Variable Models. en. Online course[END_REF][START_REF] Francesco | Predictive numerical simulations for rebuilding freestream conditions in atmospheric entry flows[END_REF].

The AM algorithm was proposed by [START_REF] Heikki | An Adaptive Metropolis Algorithm[END_REF]. The idea is to adapt the parameters of the covariance matrix of the proposal during the sampling procedure. For the first t 0 iterations, the algorithm starts as a classical MH algorithm, with a Gaussian proposal and covariance matrix K 0 . After the starting period t 0 , the covariance matrix is modified using:

K(t + 1) = K 0 , t + 1 ≤ t 0 s(M I ) K(t), t + 1 > t 0 (V.66)
The new covariance matrix K(t) is computed with:

K(t) = 1 t-1 t i=1 (x (i) I - x(t) I )(x (i) I - x(t) I ) T . s(M I
) is a tuning parameter depending of the number of inputs considered (size of M I = size of x + size of x ε ).

x(t) I is the statistical mean of the t first candidates.

V.6.3.3 Convergence assessment

Regardless of the MCMC algorithm used, the generated chain will converge towards the posterior distribution for an infinite number of iterations. However, in practice, we run many steps and assess the convergence with some empirical criteria.

First, the generated chain will visually present a distribution that doesn't present any variations after some steps. This visual criterion can be used to stop the algorithm and accept the shape of the current posterior. Quantified criteria such as acceptance rate (ratio of accepted candidates over the whole chain length) and Gelman-Rubin diagnosis [START_REF] Andrew | Inference from iterative simulation using multiple sequence[END_REF] can be used to strengthen the stopping decision.

Finally, a good practice is to consider a burn-in period T b . The first T b iterations of the chain are discarded. Indeed, as the previous algorithms start from a seed x (0) I , this method allows neglecting the influence of the beginning of the chain and a possibly miss-tuned proposal distribution in the case of the AM algorithm.

V.6.4 Use MCMC algorithms with a surrogate model

The MCMC algorithms rely on the generation of thousands of candidate samples x (c) I and the computation of acceptance probabilities. These computations require the evaluation of the model M, through the likelihood function involved in the acceptance probability. When working with expensive computer models, as in CFD, it is intractable to run MCMC algorithms. Then surrogate models, as the Kriging models (section V.4), can be used to replace the true model and generate the Markov Chains. For a comprehensive example for such an approach, the reader can refer, for instance, to [START_REF] Youssef | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF]. Before doing this, one has to ensure that the error generated by the approximation of the true model with M K can be neglected in front of the error ε.

V.7 Chapter conclusion

In this chapter, the goals and issues of uncertainty quantification approaches have been covered, including the required notions on probability theory. Then, the steps of a typical UQ workflow have been described. The methods have been described so the reader can know the necessary implementations to perform such studies. First, we covered the sampling methods and the forward propagation of uncertainties using this technique. This section stresses the need for surrogate models in UQ problems dealing with complex numerical models. The methodology to construct a Kriging model has been comprehensively reviewed. The global sensitivity analysis method using ANOVA decomposition was detailed. Finally, the goals of the Bayesian inverse problem were described, allowing us to understand the difference between deterministic and Bayesian optimizations. We finally explained the MCMC methods to address Bayesian inversion in reallife problems numerically.

Overall, this chapter depicted in a general and practical way the mathematical tools and methods to perform the uncertainty quantification analysis on the immersion cooling problems of the following chapters.

VI.1 Introduction

This chapter applies the low fidelity solver ICExo to an uncertainty quantification approach.

As we've seen in the chapter III, the model is performing the simulation of both electrical and thermal features of an immersion cooling configuration. In this previous chapter, the model was described, and its performances regarding the reproduction of original experimental data was illustrated. Briefly, the overall capabilities of the model comprise the simulation of immersion cooling BTMS by solving thermal and electrical equations to estimate the temperature, voltage and state of charge evolution of the Li-ion cells under given discharging or charging conditions. The model's capabilities were assessed by comparing its results to an original experimental test case.

VI.1.1 Reminder of the experimental setup and conditions

Before setting the mathematical framework and notations of the UQ approach for this study, let's recall the different conditions of the experimental test case reproduced with the numerical model described earlier in the chapter III. The experimental set is a pack of 8 cells immersed in a coolant flowing through an enclosure. Unsteady electrical currents are submitted to the cell under two different cycles conditions: the Datasheet cycle and the Race cycle.

The experiment aims to monitor the evolution of the temperature and voltage of the cell.

The numerical results of this section focus on the temperature prediction following those two cycles. For the reading convenience, the table describing the two cycles is shown again here:

VI.1.2 Simulation under uncertainties

The contribution of this chapter is to consider the uncertainties of the immersion cooling system and assess the predictive character of the numerical model accordingly. Indeed, the immersion cooling BTMS are multi-physics systems involving various uncertain parameters.

From the description of the physical features and parameters characterizing the electrical and thermal behavior of Li-ion batteries, which was given in chapter II, we list here the parameters considered uncertain. In the framework of this chapter, seven input parameters of the model were seen as uncertain. In the next section, more physical details are given to explain the origin of the uncertainty for these parameters.

First, thermal parameters, namely the radial and axial thermal conductivities and the specific heat of the Li-ion cells involved in the energy equation, are in practice challenging to evaluate due to the complex chemical structure of the Li-ion cells [START_REF] Drake | Measurement of anisotropic thermophysical properties of cylindrical Li-ion cells[END_REF][START_REF] Wei | A comprehensive study on thermal conductivity of the lithium ion battery[END_REF][START_REF] Chang | Thermal and solid electrolyte interphase characterization of lithium-ion battery[END_REF]. Furthermore, electrical parameters such as the internal resistance and state of health of the cell feature also complex behaviors. It is commonly agreed in the literature that the heat generated by the Li-ion cells is, for the most part, due to the Joule effect, coming from the internal resistance and the electric current going through the cells [START_REF] Park | Dynamic Thermal Model of Li-Ion Battery for Predictive Behavior in Hybrid and Fuel Cell Vehicles[END_REF][START_REF] Sato | Thermal behavior analysis of lithium-ion batteries for electric and hybrid vehicles[END_REF]. Consequently, the internal resistance parameter is expected to significantly influence the temperature evolution of the Li-ion cells when solving the heat equation. In this context, the contribution of this chapter is to use this numerical model to reduce the computational constraints inherent to CFD models and illustrate the interest in a UQ-driven approach to assess the model predictivity. Including uncertainties coming from the Li-ion batteries' thermal and electrical parameters permits assessing the impact of those multi-physical parameters uncertainties on the temperature predicted by the numerical model. Experimental data are directly used to reduce these uncertainties and learn about the parameter values and the model behaviour thanks to the resolution of a Bayesian inverse problem. Finally, a sensitivity analysis is performed to get further insights into the whole system's behaviour.

VI.2 Treatment of uncertainties

In this section, firstly, we give some elements about the physical reasons that lead to consider some parameters as uncertain. Then, having covered the uncertainty quantification mathematical framework in the last chapter, we describe these tools for the UQ analysis for the scope of this chapter.

VI.2.1 Source of uncertainties

As described in previous sections, solving the heat transfer in an immersion cooling configuration is a multi-physics problem involving physical parameters of different natures. Unfortunately, some of these parameters are hard to assess with good accuracy and confidence. Here, we propose an Uncertainty Quantification approach to address this problem. This section is devoted to characterizing each uncertainty, i.e. considering a physically sound range of variation for each parameter and associating a prior distribution. Let's review which parameters are considered uncertain in this specific case of the Li-ion batteries immersive heat transfer problem.

First, the thermal parameters, λ r , λ z , C p , are uncertain in the heat equation in the Li-ion cells, namely the thermal conductivities and the specific heat. The Li-ion cell is a "jelly-roll" layered structure of positive and negative electrodes, separator and electrolyte. A Li-ion battery generally comprises a positive electrode made of lithium cobalt oxide, manganese oxide, or nickel oxide. The negative electrode is made of hard carbon or graphite. An electrolyte of a non-aqueous solvent is intercalated between the electrodes. This structure is responsible for the anisotropic thermal behavior of the batteries [START_REF] Huat | Computational fluid dynamic and thermal analysis of Lithiumion battery pack with air cooling[END_REF]. Due to the complex chemical structure of Li-ion cells and electrochemical reactions occurring between the "jelly-rolled" layers composing the cells, the radial and axial thermal conductivities and the specific heat capacity parameters are dependent on the electrochemical state of these layers [START_REF] Drake | Measurement of anisotropic thermophysical properties of cylindrical Li-ion cells[END_REF][START_REF] Wei | A comprehensive study on thermal conductivity of the lithium ion battery[END_REF]. The highest fidelity models used to evaluate the axial and radial thermal conductivities of the batteries solve the electrochemical reactions within the layers of the cell, seen as porous media [START_REF] Wei | A comprehensive study on thermal conductivity of the lithium ion battery[END_REF].

Even using this approach, it remains challenging to assess the thermal conductivity of all the layers in the structure.

Also, the heat generated by the Li-ion cells is mainly due to the Joule effect, coming from the electric current going through the cells and its inherent internal resistance [START_REF] Sato | Thermal behavior analysis of lithium-ion batteries for electric and hybrid vehicles[END_REF][START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF][START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF][START_REF] Park | Dynamic Thermal Model of Li-Ion Battery for Predictive Behavior in Hybrid and Fuel Cell Vehicles[END_REF]. The above-described complex structure of a Li-ion cell is responsible for the obstruction of electrons transfer between the anode and cathode materials of the Li-ion cells.

The behavior of the internal resistance parameter R 0 plays a crucial role in the healing process of the Li-ion cells. Then, to represent the heat transfer in a BTMS, it is essential to obtain an accurate model or value for the resistance parameter. The internal resistance depends mainly on the temperature of the Li-ion battery and its state of charge (SOC). In practice, it is costly to assess the effects of the temperature and the state of charge on the resistance, especially for an extensive range of temperatures. Most of the models include dependency either on temperature [START_REF] Park | Dynamic Thermal Model of Li-Ion Battery for Predictive Behavior in Hybrid and Fuel Cell Vehicles[END_REF][START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF][START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF] or on both temperature and state of charge, but often for sparse values of these parameters [START_REF] Karimi | Thermal management of lithium-ion batteries for electric vehicles[END_REF][START_REF] Li | Effects of temperature on dynamic characteristics of li-ion batteries in electric vehicle applications[END_REF][START_REF] Wang | Calculation methods of heat produced by a lithium-ion battery under charging-discharging condition[END_REF]. Then it is practically challenging to assess a value for the internal resistance and its behavior in a wide range of temperatures with reasonable accuracy.

The battery's state of charge quantifies using a unitless number the amount of electric current that the battery can still deliver during a given time. It is computed from its nominal capacity expressed in [A•h] and the capacity at a given time during a charging or discharging sequence. The state of charge is a percentage of the state of health (SOH) of the battery, which is the actual value of its capacity when considered fully charged. The state of health is a way of quantifying the age of the battery. Many parameters could influence its values for a given battery, and various research efforts are focused on forecasting the evolution of this parameter [START_REF] Li | State of health estimation for Li-Ion battery using incremental capacity analysis and Gaussian process regression[END_REF][START_REF] Lu | Li-ion battery capacity cycling fading dynamics cognition: A stochastic approach[END_REF]. Considering the simulation of a battery pack with numerous cells, it seems a reasonable approach to admitting a range of variation for the state of health, in percentages, to track its influence when modeling a BTMS.

Finally, the model representing the BTMS used in this case simulates the heat transfer between the Li-ion cell and the cooling fluid flowing around. The convective heat transfer performances between the solid and fluid domains, in this case, depends on the mass flow rate ṁ and heat transfer coefficient of the fluid h f . Those parameters are also hard to assess with reasonable accuracy in practice. The coolant used in this application is a complex fluid presenting excellent heat transfer performances. However, this fluid's heat transfer co-efficient depends on temperature and pressure during operation and the flow speed used in the experiment, in a non-linear relationship. Besides, the value for the mass flow rate used during the operation of the BTMS is given indirectly through a pump generating the flow, which drives the forced convective heat transfer. Then, the value of the mass flow rate is never known precisely.

Overall, the seven following parameters are considered as uncertain: the mass flow rate ṁ, the heat transfer coefficient of the coolant h f , the internal resistance R 0 , the state of health of the cell SOH, the radial and axial thermal conductivities of the Li-ion cell λ r ,λ z , the specific heat capacity C p .

VI.2.2 Uncertainty quantification problem

Let us define the mathematical framework allowing us to deal with the uncertainties previously described. The uncertain inputs are stored in a vector X = [X 1 , ..., X 7 ]. To represent the uncertainties, each input parameter X i is seen as a random variable, taking its values according to a given distribution π i . From the physical description of the uncertainty given above, there is no prior information on the behavior of the parameter. So we consider uniform distributions for all of the seven parameters, with a range of variations detailed in Tab. VI.2, based on the user experience and the references mentioned above. The full numerical model, described in the chapter III, is from now on seen as a black box function M, giving the output from its evaluation at some input X. Here, the purpose of the study is to estimate the cell's temperature at given times of the simulation. The output of this black box computational model M is a vector of the temperatures of interest T qoi (t) evaluated at several times. Let us define the output vector of p elements as T = [T 1 , ..., T p ] with T i the temperatures at the times t i , such that T i = T qoi (t i ). The numerical model evaluation on the input X results in the output vector T according to

X

T = M(X). (VI.1)
One Uncertainty Quantification objective aims to compute some statistics on the outputs of interest. The current study sets the attention on the first and second-order statistical moments of the output vector components, namely the mean and variance. For each component T i of the vector of interest T those quantities are denoted as:

Statistical mean: E[T i ] Statistical variance: V[T i ] (VI.2)
Practically, the Monte Carlo method is used to compute those statistics on the outputs of interest. From a large number N M C of input samples X (k) , k ∈ [1, N M C ], the numerical model is evaluated through samples of the output:

T (k) = M(X (k) ) = [T (k) 1 , ..., T (k) p ] (VI.3)
Then, the statistical mean and variance of each output T i are evaluated using:

E[T i ] = 1 N M C N M C k=1 T (k) i V[T i ] = 1 N M C N M C k=1 T (k) i -E[T i ] 2 (VI.4)
To perform these computations, each input variable X i is sampled according to a uniform distribution π[X i ] described in Tab. VI.2:

X i ∼ π[X i ] = U X min i , X max i , i = 1, ..., 7
(VI.5)

VI.2.3 Surrogate model construction

To reach a proper convergence of these statistics, the number of samples N M C should be high enough. The convergence requires numerous evaluations of the computational model, which turns out to be costly. The way around consists in building a mathematical function M K , mapping the output T with the input X:

M K : R 7 → R p X → T (VI.6)
Here, we build this function using a Kriging surrogate model. To construct the model, we build a set of N LHS points called the Design of Experiments (DOE). The Design of Experiment is obtained by sampling N LHS vectors of inputs and by running the numerical model for each of these inputs. The Design of Experiments is the set of points denoted as

T DOE = X (i) , T (i) , i = 1, ..., N LHS (VI.7)
To generate the Design of Experiments, the inputs points are sampled using the Latin Hypercube Sampling technique [START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF]. Here the Design of Experiments is made of N LHS = 170 samples, chosen after a convergence study on a quantity of interest. Regarding this convergence study, the quantity of interest to assess the convergence of the design of experiment is the maximal temperature taken from the p values of the output vector. Two Design of Experiments are built constructed using respectively 170 and 420 samples. Then two surrogate models are constructed from these Design of Experiments, and this scalar quantity of interest is propagated using both models. We assess the convergence by examining the resulting

VI.2.4 Bayesian calibration problem

One problem when doing Uncertainty Quantification is the formulation of appropriate prior distributions of the input parameters. The issue is to obtain more informative distributions for the input parameters, which give a model response close to the experimental data and with reduced variability in the output predictions.

The experimental measurements of temperatures T exp can directly be used to learn more informative distributions of the inputs by solving a Bayesian Inverse problem [START_REF] Kennedy | Bayesian calibration of computer models[END_REF]. Following the Bayes paradigm, the objective is to compute the distributions of the inputs by taking into account the information provided by the experimental data. It is called the posterior distribution, and it can be computed using the Bayes theorem as seen in the section V.6.

Taking the notations of the current chapter, the Bayes formulation reads:

π[X|T exp ] = π[T exp |X] • π[X] π[T exp ] (VI.

VI.3.3 Forward propagation of calibrated distributions on a simulated race cycle and sensitivity analysis

This section aims to use distributions and associated probabilities for the input parameters from the previously described experimental conditions in a pure design perspective. Here, the input conditions for the electric current simulate an actual race cycle, composed of numerous alternating charge and discharge sequences. These sequences represent a real usage of the vehicle submitted to racing conditions. The simulated time is 1115 [sec].

The discharging sequences stand for the short and numerous accelerations of the car.

The charging sequences represent the regenerative braking that loads the battery pack when the car slows down. The alternation of these sequences responds to some Exoes' customer needs. The customer, a car manufacturer, plans to operate a racing car on the circuit for a given duration and desires to assess if the BTMS can maintain the batteries in safe ranges of temperatures during the whole racing program. Note that an additional pre-processing part is added in the numerical model in this perspective. The electrical current submitted to the batteries is computed directly from the power instructions of the electrical engine, fulfilling the manufacturer specifications for the race program. This part is off the scope of the current study. The contribution here starts from the resulting electrical current values, depending on time, which is inputted in the model as previously described in the manuscript.

The electrical part of the numerical model is adapted to represent the voltage and state of charge of a full battery pack composed of thousands of Li-ion batteries. The equations used are Eq. III.13 and Eq. III.11. The cells' connections within the battery pack allow computing the equivalent voltage and electric current, using Kirchoff's laws. For the thermal part, the heat transfer is still computed for an equivalent single-cell Li-ion representing the battery pack as described in the section III.2.1 (chapter III).

A surrogate model is built to map the inputs X to 56 temperatures of interest covering the simulated time. We use the same procedure as in sections VI. denoted:

Y 1 = min t∈[t 0 ,t f ] SOC(t) Y 2 = min t∈[t 0 ,t f ] V (t) (VI.12) Y 3 = max t∈[t 0 ,t f ] T qoi (t)
Practically, another surrogate model is constructed to map the inputs X to the three scalar indices for all three quantities of interest. Then, input parameters do not feature a strong correlation effect.

About the minimum of SOC, the parameter showing the most significant impact on its variability is the battery's state of health with a total and first-order Sobol index value close to one. For the minimum voltage, the state of health has the most significant influence on its variability. The battery's internal resistance has a non-negligible impact as well, with Sobol' indices around 0.13. Physically, the state of charge of the Li-ion cell is directly related to its voltage, so it was expected to see the SOH having the most considerable influence on this parameter. In the same way, the internal resistance is also physically directly related to the voltage. However, the Sobol analysis shows that the impact of the resistance is smaller than the state of health on the voltage variability. the internal resistance R 0 , the Joule effect seems to be the critical element governing the thermal behavior of the battery.

VI.4 Chapter conclusion

In this chapter, based on the model and the experimental test case described in chapter III, a simulation under uncertainties of the immersion cooling system was performed. Uncertainty The Sobol sensitivity analysis showed the role of each uncertain input in the variability of the several quantities of interest considered. First, the total and first-order Sobol indices showed very similar values for the three quantities of interest. Then the uncertain parameters considered don't show a strong dependency between each other in the variability of the outputs. Also, the sensitivity analysis showed that the SOH is the most influential parameter for the variability of the minimum voltage and SOC with values around 0.9 and 0.8, respectively. For the maximum temperature, the internal resistance R 0 is the most important parameter with a total index close to 0.6. The mass flow rate and heat transfer coefficient show a significant influence with total indices around 0.2.

This sensitivity analysis gives further insights related to a battery thermal management system simulation. For the thermal part, it is crucial to have an accurate model to represent internal resistance because of its impact on the battery temperature evolution. Also, representing the conjugate heat transfer by developing 2D or 3D CFD models for the fluid domain could provide more information about assessing the value of the heat transfer coefficient and guess the effect of the mass flow rate. Finally, for the electrical parameters, the state of health parameter should be simulated with reasonable accuracy to correctly represent the effects induced by the ageing of the batteries. As shown in this chapter, the SOC of the battery depends strongly on its SOH. And some references as [START_REF] Li | Effects of temperature on dynamic characteristics of li-ion batteries in electric vehicle applications[END_REF][START_REF] Chen | An internal resistance estimation method of lithium-ion batteries with constant current tests considering thermal effect[END_REF] claims the impact of the SOC on the internal resistance parameter. Developing enhanced models for these parameters should allow representing the thermal and electrical behavior of Li-ion batteries with higher fidelity.

VII.1 Introduction

From the last chapter, we have shown that the transient temperature prediction was strongly dependent on the internal resistance parameter through the sensitivity analysis using the low fidelity model. The internal resistance plays an important role in the source term of the heat equation within the batteries. This chapter aims to concentrate the efforts on a comprehensive description of this parameter. The CFD model is used here to accurately predict the immersion cooling configuration temperature. A model for the internal resistance is constructed following physical constraints and literature references. In this chapter, the test case presented in the section IV.4 from chapter IV is used to perform this analysis on the resistance model. Specifically, the test case measures the temperature at the surface of batteries immersed in air, heated by a constant electrical current.

Using an uncertainty quantification approach, the main objective of this chapter is to represent with reasonable precision the uncertainties inherent to the internal resistance parameters and evaluate their impact on temperature prediction. Also, the experimental temperature measurements provided in the present test case will reduce the uncertainties on the resistance model parameters and the resulting temperature prediction by solving the Bayesian calibration problem.

From a more practical point of view, the methodology proposed here allows reconstructing the shape of the internal resistance model indirectly from measurements of temperatures, with confidence intervals resulting from the uncertainties on the parameters. This approach might be of interest for battery systems designers as temperature measurements in battery packs are more feasible in practice than direct measurements of resistance on the batteries.

Finally, going back to a deterministic framework, we enrich the internal resistance model by considering its dependency on the state of charge. A model is constructed directly from resistance measurements at different SOC and temperature values by interpolating the data using Kriging. However, this simple approach allows assessing with small implementation effort the impact of the SOC on the resulting temperature of the heated batteries through the CFD solver. the random vector giving the Bezier construction points of the R(T ) curve. We will discuss later the choice for the distributions associated with the random vector R and then to each of the random variables R i .

VII.3 Uncertainty forward propagation

Now that we can represent the uncertainties from the internal resistance model, it is time to assess their impact on the transient temperature prediction. As described in previous chapters of the manuscript, we perform an uncertainty propagation with a Kriging surrogate model to replace the expensive CFD simulations. The definition of the model and its validation are covered in the following. The Kriging model is constructed from a Latin Hypercube Sampling design of experiments. The LHS method was used to sample N LHS = 153 models for the internal resistance

VII.3.1 Definition of the surrogate model

R (i) (T ), i ∈ [1, N LHS ].
Then, a CFD simulation is performed for each input R (i) , providing the output vector of temperature values T (i) . The points obtained after this process, constituting the Design of Experiments T DOE are denoted as:

T DOE = R (i) , T (i) , i =

VII.3.2 Choice of the prior distributions

Excepted the physical constraints detailed above on the shape of the resistance model, no a priori information is known on the behavior of the resistance model. For this reason, the introduced random variables R i are assumed to follow a uniform distribution. A physically sound interval of variation is considered from which the resistance values can be obtained by sampling the random variables following the uniform distribution: The amplitude of these intervals was determined from several trial runs on the CFD simulation, with different expressions of R(T ) models. As described earlier in the chapter IV in the test case description, the L2 error between the experimental measurements and the predicted temperature was assessed for those different trial runs. The source term expression leading to a sufficient small L2 error was selected. This so-called candidate polynomial is denoted R candL2 (T ). Supports of uniform distributions were chosen such that the CFD runs can overlap with the experimental data at all times of the simulation. More precisely, this approach led to selecting a 30% amplitude around the candidate polynomial for each of the four Bezier points. So, we have then fully defined the prior distributions for the input random vector R = (R 0 , R 1 , R 2 , R 3 ).

VII.3.3 Results of forward propagation

From these prior distributions, we sample using Monte Carlo large numbers of realizations of the random vector R. For each of these samples, a run with the surrogate model M K is performed. Then, we compute the mean and variance of the temperatures as described in The CFD run with the polynomial R candL2 is also overlapped.

This forward propagation result shows that our choice of prior distributions covers the experimental measurements and error envelopes for the full simulated time. Also, the mean of the surrogate response from the prior distribution overlaps well with the CFD model response. Then, we can assess that this modeling approach is quite satisfying regarding the ability of the model to reproduce the experimental data. The prior uncertainties assumed for the internal resistance seem consistent with the problem's physics and the CFD response.

Here we highlighted the interest of taking into account the input uncertainties from the resistance, compared to a deterministic approach using the CFD model and a single input resistance model. Indeed, the surrogate model response variability is non-negligible when the prior uncertainties on the internal resistance model are considered. Furthermore, if the mean response of the surrogate model seems to be in fair agreement with the experimental measurements, the response variability increases with time as the temperature reaches higher values. Indeed, one can note that the confidence interval amplitude increases with To build this surrogate model, the considered Design of Experiments is denoted as:

T E L2 DOE = R (i) , E (

VII.4.2 Bayesian calibration

The second approach is to consider a Bayesian calibration of the parameters [START_REF] Kennedy | Bayesian calibration of computer models[END_REF]. The objective of the Bayesian calibration is to compute the posterior distribution. In this case, the prior distribution of the input parameters was defined in section VII.2.3. There is no a priori information on the input values, except an acceptable range of variation, and then a non-informative prior (uniform distribution) was chosen. The posterior distribution is computed using a Markov Chain Monte Carlo algorithm, the Adaptive Metropolis, as described in section V.6.3.

For computational costs reasons when running the MCMC algorithm, the experimental data selected for the calibration is constituted of only 12 points. The 12 points also cover the full simulated time from t 0 to t f , but with sparser data. Then a specific surrogate similar to respectively. These reductions are computed as the ratio of the variance of the prior and posterior distributions of the inputs.

M K is built,
Also, the outcomes of the deterministic calibration are illustrated in this figure. The difference between the two methods is visible here. The L2 error minimization process gives deterministic values of resistance for the four points (green squares). Running the CFD simulation with such inputs should lead to a temperature prediction close to the experimental data, in the sense of the L2 error we have defined. On the other hand, the Bayesian calibration gives the distributions of the input values. Rather than selecting a deterministic values for each best R i , one can compute probabilities using these posterior distributions. We can then affirm that the values of the input resistance parameters will occur in a certain range, with a probability that we are now able to assess. In that sense, the Bayesian calibration is a more robust approach.

Furthermore, the Bayesian inference allows gaining knowledge on the numerical model re- ) is reduced by 98% using the posterior distributions of input parameters. Also, one can note that the variance reduction gets more significant as time in-creases. Indeed, the predicted variance computed from the prior increases significantly as the temperature increases. On the other hand, the calibrated posterior distributions allow maintaining a very small rise of the predicted variance through time. Thus, even if the variance computed from the posteriors increases slightly, the variance reduction increases as time and temperature rise. From a design point of view, this behavior is interesting since the gain of confidence in the temperature prediction increases over time. First, one can notice that the mean of the model response is slightly different from the prior predictions. The calibration changed then the bias of the prediction by the model. The average response of the model was then deviated due to the posterior distributions learned from the experimental data. Furthermore, the 95% confidence interval plots highlight the uncertainty reduction outcomes stated earlier with the distribution plots. Indeed, the interval still gets wider as the time increases, but in a significantly reduced fashion than the predictions with the priors. We can assess the temperature prediction with reasonable accuracy from this calibration process as the uncertainty is maintained in a small range for the whole simulation time. Because of the parametrization with polynomial P 53 , the temperature evolution follows the trend of the curve with R L2 (T ). Indeed the resistance field is 'flat' when the SOC is still above 50%, and the effect of SOC is not significant in this zone (Fig. Physically, we know from data that the battery's internal resistance is degrading (significant increase) when the SOC is low. Thus, the model R(SOC, T ) allows assessing the impact of this degradation on the temperature evolution, which was not observable with a model depending on the temperature only.

VII.6 Chapter conclusion

In this chapter, we illustrated a method for developing an internal resistance model, mainly responsible for heat generation. A methodology based on Bezier parameterization was proposed to build an internal resistance model dependent on the temperature only, fulfilling physical constraints acknowledged in the literature.

An accurate and fast to compute surrogate model was built to consider the uncertainties associated with the parameters driving the internal resistance model. A Bayesian inverse problem is solved using the experimental measurements of temperature directly. Informative and narrow posterior distributions are obtained for the uncertain resistance input parameters.

Thus, the uncertainty coming from these parameters is considerably reduced compared to the prior distribution assumed a priori. For the input parameters, this process allows a maximum of 76% uncertainty reduction. These posterior distributions are propagated through the numerical model, using the surrogate model. With these informative distributions, uncertainty on temperature prediction is reduced by at least 95% all along the simulated time. Overall, we have illustrated a method to assess the behavior of an internal resistance model with reasonable accuracy, indirectly from measurements of temperatures of the heated batteries. Then this accuracy was also quantified as we could assess the according to uncertainty in the temperature prediction. Enhancing the predictive character of this approach would require the use of the Bayesian framework, and the methodology presented here should be used to calibrate the parameters governing the bi-variate polynomial resistance model. Additional experimental measurements of batteries' temperature in full discharge cases would also be required to obtain the distributions of parameters representative of the overheating induced by the state of charge.

sense, it represents a valuable simulation tool for one who seeks to design an immersion cooling system and assess its performances from an overall perspective. Also, the comparison with experimental data from literature showed the good capabilities of the code to solve this kind of conjugate heat transfer problem. The analysis of the results compared to the experimental measurements highlighted the importance of the source term selection and calibration. Also, a comparison of the 2D and 3D simulations on the same test case showed which features can be represented under these two approaches.

High fidelity CFD

VIII.1.2 Uncertainty quantification methods applied to the immersion cooling of Li-ion batteries

After performing the deterministic simulations of the immersion cooling problems, uncertainties related to the physical unknowns involved in such systems have been considered. Classical uncertainty quantification methodologies, described in chapter V, have been applied to both low and high fidelity numerical solvers. We show the temperature evolution of the battery measured at mid-height on its surface, in

Overall uncertainty propagation analysis with

  Premièrement, un modèle numérique rapide, dénommé ICExo, développé chez Exoes, a été adapté pour simuler des batteries immergées dans un fluide. Le modèle numérique permet de prédire l'évolution de température et de quantités électriques (tension, état de charge) des batteries en suivant des conditions de charge et décharge variables.ixLIST OF FIGURESI.1 Share of greenhouse gas emissions by transportation type. Figure extracted from [104]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.2 Rise of patents families related to liquid immersion cooling technology over the last 70 years, presenting a compound annual growth rate (CGAR) around 40%. Data extracted and plotted from the European Patent Office database by the company Capax Infinity. . . . . . . . . . . . . . . . . . . . . . . . . . II.1 Ragone diagram including several types of batteries and combustion engines. Figure courtesy of [29]. . . . . . . . . . . . . . . . . . . . . . . . . . II.2 Jelly-roll layered structure of a cylindrical 18650 Li-ion cell. Figure courtesy of [122]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.3 Operational scheme of a Li-ion battery during discharge. Figure courtesy of [155]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.4 Relationship between OCV and SOC. Figure courtesy of [33]. . . . . . . . . II.5 Capacity of the battery against cycles. Figure courtesy of [78]. . . . . . . . II.6 Voltage evolution with SOC variation at different C-rates. Figure courtesy of [78]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.7 Geometries of Li-ion batteries. . . . . . . . . . . . . . . . . . . . . . . . . . II.8 Equivalent circuit models. Figure courtesy of [144]. . . . . . . . . . . . . . II.9 3D fields of current densities, resulting volumetric heat and temperature distributions in the case of a cylindrical Li-ion cell. Figure courtesy of [11]. II.10 Cell, module and pack levels for prismatic (a) and cylindrical (b) cells. Figure courtesy of [156]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.11 Battery pack in an immersion cooling battery thermal management system.

Figure

  Figure provided by e-Mersiv (Exoes branch). The batteries are the stacked orange cylinders. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.12 Battery module and associated cooling circuit for immersion cooling technology. Figure provided by Exoes and adapted. . . . . . . . . . . . . . . . . II.13 PCM cooling system with air channels and heat transfer principle at the cell level. Figure adapted from [118] and [83]. . . . . . . . . . . . . . . . . . . II.14 Types of indirect liquid cooling systems. . . . . . . . . . . . . . . . . . . . . II.15 Temperature field of the batteries, cooling plate and fluid in the channel. Figures from [40]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.16 Immersion cooling module developed by Exoes and e-Mersiv. . . . . . . . . II.17 Air immersion cooling device investigated in [123]. . . . . . . . . . . . . . II.18 Temperature and volume fraction of vapor in a battery pack submitted 5C discharging current. Figure courtesy of [146]. . . . . . . . . . . . . . . . . . II.19 Cooling of a hot solid by forced convection. Figure adapted from [18]. . . .

Outline I. 1 1 I. 1 . 1 1 I. 1 . 2 2 I. 2 3 I. 2 . 1 3 I. 2 . 2 4 I. 2 . 3 7 I. 3 8 I. 3 . 1 8 I. 3 . 2 8 I. 3 . 3 9 I. 4 11 I

 111111222321322423738318328339411 Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Clean mobility and challenges with Li-ion batteries for Electric Vehicles . . Battery thermal management systems . . . . . . . . . . . . . . . . . . . . . Immersion cooling of Li-ion batteries: state of the art and challenges . Challenges in battery thermal management systems . . . . . . . . . . . . . Simulation of immersion cooling systems for Li-ion batteries . . . . . . . . Uncertainty quantification for immersion cooling of Li-ion batteries . . . . Motivations and contributions of the thesis . . . . . . . . . . . . . . . . . Industrial collaboration with Exoes . . . . . . . . . . . . . . . . . . . . . . Thesis motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Scientific contributions of the thesis . . . . . . . . . . . . . . . . . . . . . . Outline of the manuscript . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure I. 1 :

 1 Figure I.1: Share of greenhouse gas emissions by transportation type. Figure extracted from [104].

Fig. I. 2 .

 2 This trend illustrates the interest of industrials in the development of this cooling system. The Li-ion batteries are directly immersed in a cooling fluid flowing in the module in this technique. The fluid goes through an entire cooling circuit, retrieving and exhausting the heat produced by the numerous batteries by conjugate heat transfer.

Figure I. 2 :

 2 Figure I.2: Rise of patents families related to liquid immersion cooling technology over the last 70 years, presenting a compound annual growth rate (CGAR) around 40%. Data extracted and plotted from the European Patent Office database by the company Capax Infinity.

  A comprehensive uncertainty quantification study was performed, relying on the predictions of the low fidelity solver. The perspective offered by this model allowed us to consider various multi-physics parameters as uncertain involved in the estimation of temperature. The characterization of uncertainties was based on literature review and the knowledge brought by Exoes experts.

Fig. II. 1 .

 1 This type of diagram plots the specific energy against the particular power for different battery systems. According to many Ragone diagrams available in the literature, Li-ion batteries stand as a promising solution to reach such objectives. This can be seen in Fig. II.1 which shows the graph of the specific energy [W • h • kg -1 ] against the specific power [W • kg -1 ] in a log scale. On the graph are plotted the capabilities of different batteries system regarding these two quantities.

  Fig. II.2 shows the typical structure of a Li-ion cell and its components for the case of a cylindrical battery. Let's review the role and a description of each of these components.

Figure II. 1 :

 1 Figure II.1: Ragone diagram including several types of batteries and combustion engines.Figure courtesy of [29].

Figure II. 2 :

 2 Figure II.2: Jelly-roll layered structure of a cylindrical 18650 Li-ion cell. Figure courtesy of [122].

  Fig. II.3. They traduce the physical process occurring in the battery [78]. The lithium cations travel in the electrolyte during discharge and intercalate into the positive insertion electrode. The positive electrode is then an electron acceptor, and the negative electrode is an electron donor. The cell voltage is a function of the concentration of the cation in the electrodes. The end of discharge is reached when the negative electrodes are emptied from Lithium cations.

Figure II. 3 :

 3 Figure II.3: Operational scheme of a Li-ion battery during discharge. Figure courtesy of [155].

Figure II. 4 :

 4 Figure II.4: Relationship between OCV and SOC. Figure courtesy of [33].

Figure II. 5 :

 5 Figure II.5: Capacity of the battery against cycles. Figure courtesy of [78].

Figure II. 6 :

 6 Figure II.6: Voltage evolution with SOC variation at different C-rates. Figure courtesy of [78].

  (a) Pouch Li-ion batteries. Figure courtesy of Flash Battery Srl. (b) Prismatic Li-ion batteries. Figure courtesy of DNK Power Co., Ltd. (c) Cylindrical Li-ion batteries. Figure courtesy of Shenzhen Himax Electronics Co., Ltd.

Figure II. 7 :

 7 Figure II.7: Geometries of Li-ion batteries.

Four main equivalent circuit

  models have been considered: the internal resistance model (Rint model), the resistance capacitor model (RC model), the Thevenin model, the PNGV (Partnership for New Generations of Vehicles) model [144, 58]. They are based on different hypothesis and then represent the behavior of the battery at a different level of fidelity as reviewed hereafter. They are illustrated in Fig. II.8: panel (a) is the Rint model, (b) the RC model, (c) the Thevenin model and (d) the PNGV model.

Figure II. 8 :

 8 Figure II.8: Equivalent circuit models. Figure courtesy of [144].

Figure II. 10 :

 10 Figure II.10: Cell, module and pack levels for prismatic (a) and cylindrical (b) cells. Figure courtesy of [156].

1 .

 1 Keep the batteries in a safe range of temperatures. The good operating range of values is commonly around 20 and 45 o C [144].

Figure II. 11 :

 11 Figure II.11: Battery pack in an immersion cooling battery thermal management system. Figure provided by e-Mersiv (Exoes branch). The batteries are the stacked orange cylinders.

  section will describe the BTMS technology and review the different existing configurations. As stated above, the batteries are set up in modules and packs, really close to each other. A typical module of Li-ion batteries is illustrated in Fig. II.11. Several modules are assembled together to compose the whole battery pack powering the vehicle, as shown in Fig. II.12. As the batteries are producing some heat during operation, it is easily understandable that disposals as in Fig. II.11 would favour hot spots formation. BTMS are systems directly integrated with the battery pack in the vehicle, and they are usually designed simultaneously. Aside from the battery pack, a circuit is integrated into the vehicle to release heat retrieved by the cooling medium. As one can see in Fig. II.12, the

Figure II. 12 :

 12 Figure II.12: Battery module and associated cooling circuit for immersion cooling technology. Figure provided by Exoes and adapted.

Fig. II. 13 .

 13 Fig. II.13. We can visualize the process: the heat transfer between the cell and the PCM is made by conduction between the two materials (Fig. II.13-c). The heat is stored in the PCM while the phase is changing. The energy accumulated in the PCM is retrieved by convection with the surrounding fluid, for instance, by forced air convection in Fig. II.13-b.

  Fig. II.14 we show several types of indirect cooling configurations.

  (a) Channels attached to the battery surface. Figure from [150]. (b) Cold plate composed of channels ensuring the flow of coolant. Figure from [40].(c) Channels meandering around the batteries. Figure from[START_REF] Cao | Thermal modeling of full-size-scale cylindrical battery pack cooled by channeled liquid flow[END_REF].

Figure II. 14 :

 14 Figure II.14: Types of indirect liquid cooling systems.

Fig. II. 15

 15 Fig. II.15 highlight this issue. Following the coolant path allows here to perfectly understand the difficulties of obtaining a good temperature homogeneity with this system. In Fig. II.15a and Fig. II.15b, the fluid is entering the channel at the inlet temperature of 25 o C . As the fluid is flowing through the channel, its temperature increases because it retrieves the heat generated by the battery cells above. Then, when it reaches the end of the circuit, its temperature is too high to ensure an excellent convective heat transfer with the channel wall and then the batteries. That's why the cells near the exit are hotter than the cells at the beginning of the circuit.Also, in this case, the heat transfer is performed at the bottom of the cells. Thus, a significant gradient of temperature is visible in every cell of the pack, as shown in Fig. II.15c. The batteries are producing heat in the whole volume. However, only the bottom part is properly cooled in this configuration. The axial thermal conductivity of the cell is not sufficient to homogenize the temperature in the whole cell, and we observe this temperature rise from the bottom to the top. Other configurations exist and aim at limiting this problem. However, if the cell is cooled at 'discrete' positions around their body only, a suitable temperature uniformity will necessarily be difficult to obtain.

Figure II. 15 :

 15 Figure II.15: Temperature field of the batteries, cooling plate and fluid in the channel. Figures from [40].

  Figure II.15: Temperature field of the batteries, cooling plate and fluid in the channel. Figures from [40].

Figure II. 16 :

 16 Figure II.16: Immersion cooling module developed by Exoes and e-Mersiv.

Figure II. 17 :

 17 Figure II.17: Air immersion cooling device investigated in[START_REF] Huat | Computational fluid dynamic and thermal analysis of Lithiumion battery pack with air cooling[END_REF].

  Fig. II.11 and Fig. II.16, cooled with fluids presenting good thermal performances. Different battery modules have been developed with various geometries and options for the flow path.

  Innovative and recent technology has received a growing interest to improve the thermal performances of immersion cooling. This technique is the two-phase (boiling) liquid immersion cooling.The configuration is the same as the liquid immersion cooling. However, the difference stands in the fluid utilized. The goal is to setup the temperature and pressure conditions in the module such that the fluid is near its boiling point. Thus, in the zone of hot spots, either due to inadequate flow circulation or overheating a battery, the temperature increases, and the fluid evaporates. As vaporization is an endothermic reaction, the heat transfer is increased significantly, and the fluid retrieves a lot of the heat released by the batteries in a two-phase state. Like in the one phase immersion cooling technology, an external circuit is added to ensure the flow of the coolant. A condensation sink can be added with the heat exchangers to make sure the fluid is re-entering the module in a pure liquid form[START_REF] Al-Zareer | A thermal performance management system for lithium-ion battery packs[END_REF]. Some authors studied numerically, and experimentally the thermal performances of this technique[START_REF] Wang | Thermal performance predictions for an HFE-7000 direct flow boiling cooled battery thermal management system for electric vehicles[END_REF]. The general conclusion is that the two-phase cooling ensures a good temperature uniformity inside the battery pack, even when discharging at high C-rates. However, performances regarding the maximal temperature reached during operation are highly affected by the fraction of vapor in the module. This fact is illustrated by a CFD simulation which results are visible in Fig.II.18. The temperature (left) and volume fraction of vapor (right) fields in a battery pack are plotted for a case of discharging at 5C. Comparing these two fields highlights the high-temperature spots formation where vapour volume fraction gets close to 0.6 (a value of 1 corresponds to saturated vapor and 0 to pure liquid). Then, such systems are challenging to manage in practice as the volume fraction of each phase needs to be carefully controlled during the charging/discharging cycles of the pack.

Figure II. 18 :

 18 Figure II.18: Temperature and volume fraction of vapor in a battery pack submitted 5C discharging current. Figure courtesy of [146].

  ∞ the temperature of the fluid, following the notations of Fig. II.19.But there are some limitations to using this approach. The heat transfer coefficient h is not trivial to assess in practice, and this calculation considers the temperature evaluated in only some points in the domains.

Figure II. 19 :

 19 Figure II.19: Cooling of a hot solid by forced convection. Figure adapted from [18].

Figure II. 20 :

 20 Figure II.20: Scheme of the conjugate convective heat transfer process in two dimensions.

  P r = molec. diffusivity of momentum molec. diffusivity of heat = ηC p λ (II.16) η is the viscosity of the fluid in [N • s • m -2 ]. C p the specific heat capacity and λ the thermal conductivity, both defined earlier in the manuscript. Very small values of the Prandtl number (P r 1) indicates that the thermal diffusivity dominates over the momentum diffusivity. The energy transfer in the fluid occurs mostly due to heat dissipation controlled by the thermal conductivity. On the opposite, large values (P r 1) show that the fluid dissipates its energy mainly by the momentum dissipation, caused by a significant viscosity value.

  Fig. III.1, the general process of the model is illustrated. The electrical equations representing the ECM are coupled to the heat equation in the battery. Heat transfer between the solid and the fluid domains is computed through the balance of specific enthalpy in a conduct.

  configuration. For the thermal part, the variables computed by the model are the Li-ion battery's temperature and the surrounding fluid's temperature. The voltage and the state of charge of the battery are computed in the electrical part. The description of the whole numerical model has been summed up in Alg. 1.

(

  Figure III.1: Scheme of the ICExo model

Figure III. 2 :

 2 Figure III.2: Mesh of solid domain with temperature unknowns in the center of each mesh cell

  The equation system formed by Eq. (III.1), Eq. (III.4), Eq. (III.8) and the boundary conditions allows to compute the three unknowns of the problem: the temperature in the solid domain T (x s , t), the specific enthalpy of the fluid h(x f , t) and the temperature in the fluid domain T (x f , t), at any time t. The model's specific output of interest is the mean of temperature within the battery on the right side, near the solid-fluid interface, as shown in Fig. III.3. This zone is denoted Ω right int and its surface S Ω right int

Figure III. 3 :

 3 Figure III.3: Computational domains: solid domain and fluid layer, control volume scheme for the specific enthalpy computation.

  2 [V] between two consecutive data. Following the RC circuits wiring in Fig. III.4, the impedance Z of the battery is then expressed:

Figure III. 4 :

 4 Figure III.4: Adapted Thevenin equivalent model of the Li-ion battery.

  with 32 cells. The packs are immersed in the fluid within a tank or enclosure. The fluid utilized is either air or coolants with better thermal performances. Cells temperatures are usually measured by thermocouples placed on the skin of the cells. Using thermocouples on each cell allows evaluating the temperature uniformity within the packs. These setups focus on different parameters of interest, such as maximum temperature, maximum spatial temperature deviation and temperature uniformity, giving informative insights on the behavior and performances of immersion cooling systems. Let's review some results from the literature provided by such experimental investigations.

Figure III. 5 : 20 Table III. 1 :

 5201 Figure III.5: Experimental set up. Position of the temperature sensors (left), computer aided design (middle) and real set up (right).

Figure III. 6 :

 6 Figure III.6: Experimental measurements for the Race cycle.

  thermal prediction capabilities of the model are shown in Fig. III.8 and Fig. III.10. The temperature evolution computed by the numerical model is compared against the experimental data (Anode 1 temperature, see Fig. III.5) and measurement error envelope. The same values of input parameters are taken for both cases, so it is possible to compare the model's performances quantitatively in each situation. The values of the thermal input parameters used in these simulations are given in Tab. III.3. Those values were obtained after the calibration process detailed later in the manuscript (see chapter VI). Especially the range of

Health SOH 9 . 1 Table III. 3 :Figure III. 8 :

 9138 Figure III.8: Temperature computed by the numerical model for the Datasheet cycle case.

Figure III. 9 :

 9 Figure III.9: Temperature computed for two DataSheet cycles.

Figure III. 10 :

 10 Figure III.10: Temperature computed by the numerical model for the Race cycle case.

IV. 1 68 IV. 1 . 1 68 IV. 2 72 IV. 2 . 1 72 IV. 2 . 2 77 IV. 3 79 IV. 3 . 1 93 IV. 4 . 1 103 IV. 4 . 4 109 IV. 5 111 IV. 1

 1681168272217222773793193411034410951111 CFD model for conjugate heat transfer of Li-ion batteries . . . . . . . . TrioCFD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 IV.1.2 CFD for immersion cooling: technical literature review . . . . . . . . . . . Governing equations and numerical methods in TrioCFD . . . . . . . . . Description of the physical problem . . . . . . . . . . . . . . . . . . . . . . Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation of 3D forced convective heat transfer . . . . . . . . . . . . . Test case description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 IV.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 IV.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 IV.4 Simulation and validation of natural convection heat transfer . . . . . . Experimental test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 IV.4.2 2D simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 IV.4.3 3D simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Comparison 2D/3D and experimental data . . . . . . . . . . . . . . . . . . Chapter conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CFD model for conjugate heat transfer of Li-ion batteries IV.1.1 TrioCFD The second code used in this thesis is TrioCFD, seen as a high fidelity (HF) solver compared to the numerical model of the previous chapter. TrioCFD is a open source Computational Fluid Dynamics code developed at CEA since 1995 [140]. It is based on the Finite Element Method (FEM). Its applications are the simulation of incompressible and quasi-compressible flows and the heat transfer associated with nuclear applications. It is massively parallelized through the use of MPI libraries. In the framework of this thesis, the code is used for the first time in applications related to the conjugate heat transfer of Li-ion batteries in immersion cooling configurations. Nuclear applications and the immersion cooling of Li-ion batteries have a lot in common regarding the heat transfer problem definition. However, the magnitudes of velocity fields or temperature reached are significantly different in these two domains. This solver has been chosen to perform the simulations in the framework of this thesis, as it is open source and new developments would have been easier to implement if needed. Also, one development objective of this code is to represent multi-phase flows, being a technology of interest for immersion cooling systems. Addressing the simulation of immersion cooling of Li-ion batteries with this model requires a specific approach described in the following. The meshing of the solid (batteries) and fluid (surrounding gas coolant) is performed with unstructured grids. The full system of partial differential equations is solved in the solid and fluid domains. The electrical behavior of the cell is not solved directly as it is done in the model ICExo. However, the electrical features of the batteries are represented by considering laws for some relevant parameters, as the approach in section II.2.2.2 from chapter II. Then those electrical parameters are directly coupled to the set of equations for the fluid and the heat transfer. The electrical parameter considered is the internal resistance R and the electrical current I, playing a role in the heat source term of the solid domain. Thus, the full model allows linking the electrical features with the thermal response of the batteries. A general overview of the model operation as used in this thesis is schematized in Fig. IV.1.

Figure IV. 1 :

 1 Figure IV.1: High-fidelity CFD model

8

  Turbulence model added to RANS equations if needed. Otherwise a laminar flow is considered. 9 Simulations in 2D or 3D geometries.10 In 3D simulations with natural convection: Boussinesq approximation.11Perfect thermal contact between solid and fluid: continuity for temperature and heat flux at the interface.

  Similarly to IV.2.1.1, the energy conservation principle is applied to the fluid domain. The energy conservation equation with the variables related to the fluid domain reads:

  location and time. The fluctuating terms u , p stand for the deviations to this mean. Averaging the previous Navier-Stokes equations allows to get rid of the local fluctuations, and the system of equations considering the average and fluctuating variables now reads: ∇ • u = 0 (IV.13)

  For the energy equations, degrees of freedom for the temperature are also located in the center of the faces. An illustration of the elements and the unknown location in 2D is given in Fig. IV.2.

Figure IV. 2 :

 2 Figure IV.2: Element in two dimensions and location of the degrees of freedom

Fig. IV. 3 .

 3 The solver for the fluid and solid domains present different time steps, respectively δt f and δt s . Each solver runs independently in parallel for a cycle of duration ∆t, once the boundary condition at the interface has been updated. The coupling time step ∆t is global, linked to the coupled problem. The physics of the problem governs the choice of ∆t according to a trade-off between accuracy and precision. Following the notations of Fig. IV.3, at the coupling step t N = N ∆t the equality of flux is assumed along the interface Γ. This assumption allows to compute the contact temperature on the interface boundary, denoted T I (x, t N ), by integrating the flux equality expression over the interface:

  the interface boundary for both fluid and solid domain. Then the resulting temperature fields are computed in parallel until the next coupling step (N + 1)∆t.

Figure IV. 3 :

 3 Figure IV.3: Parallel coupling strategy. Figure courtesy of [112].

Figure IV. 4 :

 4 Figure IV.4: The 9 electric cells (solid domain Ω s ) and surrounding box representing the fluid domain Ω f

  3. The reader can refer to the Fig. IV.4 to match the boundary names to the geometry. The air flow is entering by the whole surface of the boundary ∂Ω inlet at a speed of 0.1 m.s -1 along the x axis. The flow

cell number 5 Figure IV. 5 :Figure IV. 6 :Figure IV. 7 :

 5567 Figure IV.5: Meshes and number of tetrahedras. View from above of the cell 5

  Fig. IV.8, we plot the evolution of the temperature at the center of each cell. The temperature is plotted for cells 2, 4, 5 and 8 (refer to Fig. IV.4 for the position of the numbered cells). Cells 4 and 5 arebehind each other with respect to the flow direction. Cell 4 presents a final temperature lower of 0.5 degree. We can then assume that cell 4 is better cooled than the 5. It can be explained that the cells upstream of the flow receive air that has not exchanged energy with any cells yet. However, by plotting the temperature of the cells 2, 5, and 8, which are in the same row perpendicular to the flow direction, we can assess a good temperature homogeneity between the cells in the cross direction of the flow.

Figure IV. 8 :

 8 Figure IV.8: Evolution of the temperature at the middle of the cells.

Figure IV. 9 :

 9 Figure IV.9: Temperature field [K] of solid domain and velocity vector field [m • s -1 ]

Fig. IV. 10 .

 10 Fig. IV.10.From the above view plot, we can see that the temperature presents a significant gradient in the flow direction for the cells in the front row. Again, the air is entering the inlet boundary at the low inlet temperature, and the front face of the cells can present a higher heat transfer. Also, the temperature gradients due to the recirculating area are visible in the profile view. Note that the top of the cells presents lower temperatures. Due to the geometry, the airflow presents high speeds at the top of the domain, and we can expect a good heat transfer in this area.To confirm the fact that the heat transfer is decreased in the low air velocity zones, we plot the two temperature fields in Fig. IV.11. Note that the temperature scales are different from solid and fluid domains to see the temperature gradients in the solid domain more precisely. The zones where the velocity is near zero in Fig. IV.10 match with the zones presenting the higher temperatures in the fluid domain. As in these areas, the temperature of the solid and fluid are closer, the convective heat transfer rate decreases and the solid domains present the hot spots. Similarly, the zones with high-speed velocity match with the low-temperature zones.

  Figure IV.11: Temperature fields [K] in the fluid and solid domains.

  Fig. IV.12, Fig. IV.13, right). In Fig. IV.12, the temperature plots along the cross flow direction are shown. The precise location of the monitored temperature is specified by the line in the geometry scheme on the right. The temperature is monitored at three different heights in the z direction. The direct comparison of FLUENT and TrioCFD is made by overlapping the temperature curves. Be aware that the temperature is plotted in Celsius degrees to stick with the reference results in these figures. First, this way of presenting the results supports the previous observations about the hot spots. From figure Fig. IV.12 one can quantify the amplitude of those hot spots in the zone behind the cells where the air velocity is relatively low. Notably, in the figure IV.12 a temperature rise from the inlet air around 3 o C can be observed in the zone behind the middle row of cells, at z = 50 mm. With a similar approach, Fig. IV.13 allows getting specific knowledge about the airflow characteristics within the fluid domain. The norm of the velocity vector is plotted in the crossflow direction. The link between the temperature hot spots and air velocity can be assessed through those different plots. By comparing the figures Fig. IV.12 and Fig. IV.13, one can see that the location of the hot spots match exactly with the location where the velocity is almost null.

Figure IV. 12 :Figure IV. 13 :

 1213 Figure IV.12: Variation of temperature at three heights in cross-flow direction, downstream (t=500 sec)

  The comparison with the experimental data should allow us to understand the heat transfer process of immersed Li-ion batteries. Unlike the previous section, the test case treated here consists of heat transfer with natural convection. In terms of flow resolution, this simple configuration should allow addressing the following issues. First, assess the accuracy of the model for the heat source of solid domain equations, supposed to represent the heating patterns of Li-ion batteries. The purely numerical comparison couldn't provide such validation in the previous test case. The experimental measurements of transient temperature evolution should validate or not the behavior of this source term model. Second, a better understanding of the heat transfer process between the fluid and the solid would be featured in this case of natural convection. The experimental test case is reproduced using 2D and 3D simulations. The comparison of these numerical predictions, coming with their respective assumptions with respect to the experimental data, may highlight the physical processes occurring in this configuration.

  experiment reproduced in this study is taken from He at al. in the International Journal of Heat and Mass Transfer [57]. The setup consists in two parallel rows of four Li-ion cells each (see Fig. IV.14). The Li-ion cells used in the experiment are cylindrical 26650 cells with LiFePO4 chemical composition. The battery pack, composed of two rows of cells, is set inside an enclosure made of plexiglass walls. Each cell present a diameter of d = 25.85 [mm] and height of L = 62.5 [mm], with a given spacing between cells illustrated in Fig. IV.14. A constant discharging electric current is applied to the Li-ion cells during the whole experiment. According to the discharge rate of 1.5C, the current across the entire battery pack is set to I = 6.9 [A]. This steady discharge current is applied to the cells from the beginning of the experiment, at time t 0 = 0 [sec] until the final time t f = 1600 [sec]. The experiment aims to monitor the temperature evolution of the heated Li-ion cells measured at key locations of each cell surface. The black cross represents the positions of the probes in Fig. IV.14b. The probes are situated at the mid-height of the batteries.

Figure IV. 14 :

 14 Figure IV.14: Experimental set up. Figure courtesy of He et al. (2014) in the International Journal of Heat and Mass Transfer, vol. 72 [57].

Figure IV. 15 :

 15 Figure IV.15: Mesh, computational domains and boundary labels for the CFD simulation

  Five meshes are constructed to assess the convergence. They are represented in Fig. IV.17, with a focus on cell 1 and the plexiglass wall. For each mesh, the corresponding number of triangles of the whole mesh is indicated. The monitored quantity for the convergence is the temperature at the solid/fluid interface of cell 1, at the position represented by the white cross in Fig. IV.15. We treat the case of natural convection with a null inlet velocity, using the CFD configuration described just above. For the mesh convergence only, the simulation is performed with a constant source term in the heat equation q g = cte. The simulation is performed for the whole experimental time of 1600 seconds. The maximal time step authorized for the implicit Euler scheme is δt = 5E -2 sec. The results are shown in Fig. IV.16. The convergence is reached for the meshes containing 2.8E + 4 triangles or more. In the following, results presented on this test case in 2D are obtained using the mesh with 3.4E + 4 triangles, in Fig. IV.17d. The simulations with this mesh took 40 minutes to run in parallel on 36 cores.

Figure IV. 16 :

 16 Figure IV.16: 2D meshes convergence (number of triangles). Temperature on the cell 1.

4 Figure IV. 17 : 97 IV. 4 . 2 . 3

 41797423 Figure IV.17: 2D meshes -number of triangles 97

Figure IV. 18 :

 18 Figure IV.18: Temperature evolution for the four probes location.

  presented, i.e. with air velocity set to u = 0 [m • s -1 ] at inlet boundary. A Dirichlet boundary condition is applied to the temperature of the inlet boundary: T inlet = 295.2 [K]. The same value is used for the initial condition in the whole domain.

  They are represented in Fig.IV.[START_REF] Ceraolo | Experimentally-Determined Models for High-Power Lithium Batteries[END_REF] for different times, using a color scale covering the range of temperatures from t 0 to t f . This figure highlights the process of transient conjugate heat transfer. The heat is generated in the solid domain Ω s with the volumetric heat source q g from the heat equation. Then the heat is spread through the solid domain. The fluid in Ω f surrounding the solid domain is retrieving the heat generated through the coupling conditions. This kind of thermal behavior is expected regarding the geometry of the experimental setup. Then, this simple configuration underlies the issues of thermal distribution in real battery packs where several Li-ion cells are arranged close to each other. The present numerical model can illustrate the difficulty of evacuating heat from the inner part of the pack, where the cells are in a confined zone.

  Figure IV.19: Temperature fields evolution, in [K]

Figure IV. 20 :

 20 Figure IV.20: Temperature field at t = 1600 sec, in [K]

Figure IV. 21 :

 21 Figure IV.21: Temperature predicted by the 2D CFD model and experimental measurements.

Figure IV. 22 :

 22 Figure IV.22: Geometry and mesh in 3D

Figure IV. 23 :IV. 4 . 3 . 2

 23432 Figure IV.23: 3D meshes convergence (number of tetrahedra). Temperature on the cell 1.

  Figure IV.24: Temperature [K] and velocity [m • s -1 ] fields

  2D and 3D simulations and with respect to the experimental measurements. The results are shown in Fig. IV.26. We overlap the last 2D CFD results and the experimental measurements with the 3D results. Both simulations were run with the same source term computed through the polynomial R CL2 (T ).

Figure IV. 26 :

 26 Figure IV.26: Temperature at the probe position for the cell 1.

  The model is applied to two different applications. The first test case is a code-to-code comparison on a forced convective heat transfer with air. The results allow us to assess the convergence of the TrioCFD code on such a configuration. Also, the mechanisms of convective heat transfer are observed thanks to this simulation. The link between the cooling fluid flow and the temperature is shown. Poor heat transfer performances occur where the flow present recirculating or low-velocity zones.

  surements are usually seen as the data presenting the highest level of fidelity regarding the physical phenomena under investigation. The simulation's required work consists of setting up a model the closest to the real world experimentations and physical processes. The quality of a model is assessed with respect to its ability to reproduce the experimental data. This is the validation step of the so-called Verification and Validation process of a computer model, illustrated by [105] in Fig. V.1. However, taking the example of simulation for Li-ion batteries, we have seen in the first chapter that the accuracy of the models would depend on the scale considered to represent the physics of the immersion cooling problems. There are many different multi-physics parameters involved in the immersion cooling problem. Those parameters are a function of the inner chemical phenomena at stake in the structure of the batteries. It is a complex challenge to give reasonable confidence values for these parameters. Choices for these parameters rely on models close as much as possible to reproduce experimental data. On the other hand, the experimental results are an outcome of real-world physics, and are intrinsically uncertain due to the tolerance of the measurement tools.

Figure V. 1 :

 1 Figure V.1: Validation process. Figure extracted from [105].

  Fig. V.2, where we prefer to focus on the uncertainty characterization of the inputs to simplify the sketch of the workflow.

  Fig. V.3 that the LHS sampling technique with the sub-intervals decomposition (dark dotted lines) can cover a larger area of the square than the MC samples.

Figure V. 3 :

 3 Figure V.3: Sampling of a random vector using LHS and MC techniques.

  Fig. V.4 is illustrated a typical random process in a one dimensional case: a collection of random variable indexed by x: {Y (ω)} x∈X defined on the same probability space. A realization Y (ω) of the stochastic process is drawn and two evaluations of this process Y (x, ω) and Y (x , ω) are visible.

Figure V. 4 :

 4 Figure V.4: Illustration of a 1D stochastic process Y , indexed by x ∈ R. Figure courtesy of [14].

  Figure V.5: Plot of GP with increasing number of realizations and PDF associated to the realizations in x = x 0 . Example adapted from[START_REF] Roelants | Gaussian processes -from scratch[END_REF].

  Now, in the interpolation perspective, let's define the Kriging predictor, i.e. the expression defining the evaluation of the Kriging model on a new point x : M K (x). The goal is to predict the value of the Kriging model at a new input x, based on a set of observations of the true model. The observations are contained in the so-called Design of Experiments (DOE).

Figure V. 6 :x i -x i θ i 2 0. 5 (V. 39 )V. 4 . 3 . 3

 6539433 Figure V.6: Gaussian correlation function and resulting GP for different θ.

  3 are used to compute those integrals. Then, the results of the Sobol analysis has to be converged with respect to the number of Monte Carlo samples. Also, for computational cost reasons, a surrogate Kriging model can replace the model to perform the computation of the Sobol index. Finally, one can note that the Sobol indices values might vary with the PDF π(x) against which are integrated the model responses. The choice of these distributions is discussed comprehensively in the next section by considering the notions related to model calibration. The Bayes theorem using these distributions is expressed as: π(x I |y exp ) = L(x I ; y exp )π(x I ) X I L(x I ; y exp )π(x I )dx I (V.55) The desired distribution π(x I |y exp ) is called the posterior distribution. It gives the information about the unknowns x I after being conditioned to the observed data y exp . L(x I ; y exp ) is the likelihood function computing the probability of having observed the measurements y exp given the realizations of inputs x I . The denominator term Z is a normalization constant: Z = X I L(x I ; y exp )π(x I )dx I It is the likelihood integrated over the prior distribution.

2 . 2 •

 22 Algorithm Metropolis-Hastings algorithm Initiate Markov chain with seed x (0) I ∈ X I . for t = 1, ..., T • Draw a candidate sample x Draw a sample v from a uniform distribution U([0, 1]). • Perform an acceptance test for the candidate: if v ≤ α(x I else (i.e. v > α(x (c) I , x (t) I )) reject the candidate and x (t+1) I = x (t) I The MH algorithm calls evaluations of the unnormalized posterior density which can be 140 given directly by the evaluation of the likelihood function and the prior: π(x I |y exp ) ∝ L(x I ; y exp )π(x I ) (V.64)

  Finally, parameters related to the fluid domain present as well some uncertainties. The mass flow rate of the cooling fluid is given indirectly by the operation of the pump generating the convective heat transfer. Due to the geometry of the arrangement of the Li-ion cells within the pack, it is hard to assert a value of the mass flow rate from the pump input instructions. Overall, uncertainties from those multi-physics parameters will generate considerable variability in the numerical model's temperature, voltage, and state of charge prediction.

  distribution of the quantity of interest. The distributions computed from both models are illustrated in Fig. VI.1. The two distributions overlap over their full support, and then we can consider that the surrogate model has converged. We stick then to the Design of Experiments containing 170 samples for the surrogate models constructed later in this chapter.

Figure VI. 1 :

 1 Figure VI.1: Convergence assessment of the surrogate models, for the scalar quantity maximal temperature.

Figure VI. 2 :

 2 Figure VI.2: Plot of the values of temperatures computed by the surrogate model (Y Kriging ) against the values from the numerical model (Y N umerical model ), evaluated on an independent validation set of input points.

  8)π[T exp |X] stands for the likelihood and π[X] is the prior distribution defined in the previous section. π[T exp ] is seen as a normalization constant called the marginal likelihood. The section VI.2. The MCMC simulation allowed us to detect which values of the inputs are more likely to give a numerical response close to the experimental data, i.e. the values close to the peak of each distribution.

Figure VI. 3 :

 3 Figure VI.3: Posterior distributions of the input parameters after the calibration.

Figure VI. 4 :

 4 Figure VI.4: Distributions of the prior and posterior predictions for the six quantities of interest. Comparison with experimental measurements and error envelopes.

Figure VI. 5 :VI. 3 . 2

 532 Figure VI.5: Posterior and prior predictions of the surrogate model mean and 95% confidence interval. Comparison with experimental measurements and error envelope.

Figure VI. 6 :

 6 Figure VI.6: Posterior predictions in the Race Cycle case, surrogate model mean and 95% CI interval. Comparison with experimental measurements.

Figure VI. 7 :

 7 Figure VI.7: Posterior and prior predictions of the surrogate model mean and 95% confidence interval, for the real racing case.

Finally, the thermal

  Figure VI.9: First order Sobol' indices for the three quantities of interest.

  quantification methods, such as Kriging surrogate models and Monte Carlo evaluations, allowed treating the uncertainties coming from the model's physical parameters. Experimental data were used directly to compute posterior distributions of the input parameters by solving the Bayesian inversion problem. Posterior learned distributions significantly reduced the uncertainty in the model prediction of temperatures. The model and calibrated distributions showed good reproducibility capability when comparing the numerical results to the experimental data under different conditions. Specifically, the support of the posterior distributions is significantly reduced with respect to the prior. The propagation of the posterior distributions showed a good ability to match the experimental data, with a very narrow 95% confidence interval in both Datasheet and Race cycles. Finally, the posterior distributions and the numerical model were used to assess the temperature response of a BTMS under realistic racing conditions.

Figure VII. 2 :

 2 Figure VII.2: Five samples of the Bezier control points and corresponding R(T ) curves (matched by colors). The red dotted lines represent the temperature range of the experiment. The four abscissa of the control points are: T 0 = 290 [K]; T 1 = 294.8 [K]; T 2 = 299.5 [K]; T 3 = 304.1 [K]

Formally, the problem 2 (VII. 3 )

 23 under consideration here is the computation of mean and variance for the quantities of interest, resulting from the uncertainties modeled by the random vector R defined just earlier. The quantity of interest outputted by the model is a vector of p temperatures, at p times of the transient simulation: t 0 ≤ t 1 ≤ ... ≤ t p ≤ t f . They are stored in a vector denoted: T:T = (T 1 , ..., T p ) (VII.2)The statistics of interest are the mean and variance of the temperatures T j at fixed timestj of the experiment. Their computation from the Monte Carlo samples leads to the following expressions for the mean E[•] and variance V[•] for each output T j : The evaluation of many samples N M C to make the statistics converge properly corresponds to as many runs of the CFD model, each corresponding to an evaluation of the computational model with an internal resistance model R (k) (T ) as input. It turns out to be prohibitive with the present CFD model, which presents a high computational cost. Then we build a surrogate model of the vector quantity of interest as a function of the input parameters, which can be used instead of evaluating the CFD model. Practically, let's define for the surrogate model the input as a vector of four sampled values of resistance R = [R 0 , ..., R 3 ]. The output is represented by the vector T of p temperature of interest values Tk , at different times of the experiment, as defined just above. To represent the full simulated time t f = 1600 sec, we choose to use p = 81 values in the vector of output temperatures.

M

  K : R 4 → R p R → T (VII.4)

1 ,

 1 ..., N LHS (VII.5) The expensive part in terms of computational cost is the construction of the Design of Experiments: each of the N LHS simulations was performed in parallel on an HPC cluster PlaFRIM, using 36 CPU cores, for a computational time of around 40 minutes. So the overall cost for the Design of Experiments construction involved 153 × 36 cores ×40 minutes. Once the surrogate model is built, we should assess that it represents the numerical model accurately. The Fig. VII.3 illustrates the evaluation of temperatures by the surrogate model for some inputs R (i) versus the one evaluated with CFD. All the points overlap close to the y = x line, showing an excellent accuracy of the surrogate model. In practice, an other validation Design of Experiments T val DOE containing N val = 82 points, was used to perform this validation.

Figure VII. 3 :

 3 Figure VII.3: Comparison of surrogate model evaluations of temperature M K with the CFD model M for given inputs R (i) . The red line represents the first bisector curve.

R

  i ∼ U([R min i , R max i ]) for i = 0, ..., 3 (VII.6)

  Eq. (VII.3). The mean and variance envelope, corresponding to the 95% confidence interval of the surrogate model response, is represented by the blue dots and envelope in Fig.VII.4.

VII. 4 . 1 2 (Figure VII. 5 :

 4125 Figure VII.5: Computation of the error between CFD and experimental data

  i) L2 , i = 1, ..., N LHS (VII.9) Then a deterministic optimization problem is solved to get the values of input R which minimizes the function f K . This outcomes in an optimized vector called R L2 defining a model R L2 (T ) which gives the CFD response the closest to the experimental measurements, in Fig. VII.6 (green). The temperature prediction resulting from the CFD computation using this input resistance model can be seen in Fig. VII.7 (green).

  The objective is to compute the distribution of the input parameters R = [R 0 , ..., R 3 ] conditioned to the experimental data T exp . Practically, this is the distribution of inputs leading to the model response the closest to the experimental data in a Bayesian sense. Following Bayes theorem, the relation between the posterior distribution, the likelihood, the marginal likelihood and the prior distribution reads:π[R|T exp ] = π[T exp |R] • π[R] π[T exp ] (VII.10)where π[R|T exp ] is the posterior distribution, π[T exp |R] the likelihood and π[R] the prior distribution. π[T exp ] is the marginal likelihood used essentially as a normalization constant.

3 Figure VII. 6 :

 36 Figure VII.6: Prior (blue) distributions and posterior (red) distributions of input parameters after Bayesian calibration. In green, resulting values of resistance after deterministic calibration.

  sponse. The Fig. VII.7, shows the temperature distributions (red) predicted by the surrogate model used in the calibration, resulting from the propagation of posterior input distributions in Fig. VII.6. The posterior predicted distributions are more narrow than the priors. The predicted variance from prior and posterior distributions, together with the associated variance reduction, are given in Tab. VII.1. For instance, the uncertainty for the predicted temperature at t = 1517 [sec] (Fig. VII.7

Figure VII. 7 :

 7 Figure VII.7: Prior (blue) and posterior (red) predicted distributions of the quantities of interest. In green, result of the CFD model with the R L2 (T ) input model.

Figure VII. 8 :+ a 30 x 3 + a 21 x 2 y + a 12 xy 2 + a 03 y 3 + a 40 x 4 + a 31 x 3 y + a 22 x 2 y 2 + a 13 xy 3 + a 50 x 5 + a 41 x 4 y+ a 32 x 3 y 2 + a 23 x 2 y 3 VII. 5 . 2

 8323423542352 Figure VII.8: Mean and 95% CI with prior and posterior distributions

  .10). After t = 1000[sec], when the SOC goes under 55%, the slope of the temperature curve is increasing. One can notice that the temperature curves corresponding to R L2 (T ) and R(SOC, T ) models show different slopes from this point. After t = 1600 [sec], when the SOC goes under 30%, the slope is increasing even more. At the end of the simulation, when the SOC approaches 10%, the curve presents a local and high increase. The resistance value is in the steep and high-value zone of the field from Fig. VII.10. Then we highlighted the modified behavior of the temperature response when including the effect of SOC.

Figure VII. 11 :

 11 Figure VII.11: CFD temperature response with R L2 (T ) and R(SOC, T ) models.

Finally, an approach

  to enhance the internal resistance model fidelity is proposed. A resistance model considering the state of charge and temperature's combined effect is constructed based on sparse experimental data. The Kriging parametrization allows assessing the relevance of the model regarding the available data while keeping a smooth and realistic shape of the bi-variate function. The model is propagated through the CFD code, and transient simulations are performed for longer physical times. The model taking into account the effect of SOC shows a significant slope increase in the temperature prediction when the Li-ion cells reach low SOC values at the end of the discharge.

  solver for transient conjugate heat transfer (chapter IV) While the low fidelity solver showed good predictive capabilities, it cannot solve some of the physical processes involved in the conjugate heat transfer problems applied on the immersion cooling of Li-ion batteries. The numerical simulation for this application was then addressed with an open-source CFD solver, TrioCFD. Based on the Finite Element Method, this code was strongly validated in the past at CEA on nuclear applications. It was used for the first time on battery cooling problems during this thesis. The code was validated against another commercial software, FLUENT, on two different test cases.

  LF solver (chapter VI)The low fidelity solver ICExo features various multi-physics modeling parameters, which are uncertain inputs. Parameters related to the thermal and electrical physics of the battery and the coolant flow have been considered uncertain. The effect of the variability of these parameters on the quantities of interest has been assessed by performing forward uncertainty propagation assisted with a surrogate model. Based on the input uncertainty, the confidence intervals and the statistical mean of the temperature evolution have been computed. These statistics on the quantities of interest have been computed for different experimental conditions. Also, the sensitivity analysis allowed drawing a hierarchy of the importance of these various parameters concerning the thermal and electrical performances of the system. For the maximal temperature reached during the simulated race of the vehicle, the internal resistance is the most influential parameter. Heat transfer coefficient and mass flow rate of cooling fluid feature also non-negligible importance on the variability of this quantity. Regarding the first electrical quantities of interest, the minimum of voltage, the battery's state of health is the most determinant parameter. The internal resistance also comes into play but with less importance. Then, the minimum state of charge for the battery is the second electrical quantity under study. For this characteristic, the battery's state of health was predominant by far compared to the others. In a nutshell, the electrical quantities of interest are influenced by During the thesis, a part of the work not covered in this manuscript has been dedicated to validating this new operator's behaviour in TrioCFD. This operator was validated on an experimental test case presented in[START_REF] Drake | Measurement of anisotropic thermophysical properties of cylindrical Li-ion cells[END_REF]. It is the case of cylindrical batteries heated radially with a steady heat source. The experimental results were reproduced with TrioCFD by solving only the heat equation in a cylindrical battery represented in a 2D geometry. In this case, as stated in the chapter II, the radial conductivity is much lower than the axial one. The code gave results that compare well with the experimental data and the analytical model provided in this reference. The results of this comparison with TrioCFD are provided in Fig. VIII.1.

Fig.

  Fig. VIII.1a, The heat source is applied on the whole surface of the battery along the radial axis e r . The temperature field resulting from the radial and constant heating is shown in Fig. VIII.1b.From these results, we can expect that the anisotropic behavior can impact the temperature prediction in the resolution of full conjugate heat transfer problem. Significantly, the temperature field in the solid domain seems strongly impacted by the anisotropy. Then, the choice of the battery arrangement (aligned, staggered, etc...) in a larger battery pack could be designed considering this temperature distribution. To obtain a better heat transfer between the batteries and the coolant, the flow path of the fluid or the position of the batteries could be adjusted by considering the influence of anisotropy on the resulting temperature field in the solid domain.

Figure VIII. 1 :

 1 Figure VIII.1: Anisotropic thermal conductivity test case

EU (Convention) -Share of transport greenhouse gas emissions

			Other Transportation
			Aviation
	13.4%	13.9%	Road transport
			Maritime
			Railways
	71.7%

II.1 Lithium-ion batteries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

	II.1.1 Generalities and definitions . . . . . . . . . . . . . . . . . . . . . . . . . . 14
	II.1.2 Electrochemical reactions . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
	II.1.3 Electrical characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
	II.1.4 Types of Li-ion batteries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
	II.1.5 Thermal issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

II.2 Modeling approaches for Li-ion batteries and heat source . . . . . . . . 24

	II.2.1 Modeling the characteristics of Li-ion batteries . . . . . . . . . . . . . . . . 24
	II.2.2 Heat source term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

II.3 Immersion cooling of Li-ion batteries . . . . . . . . . . . . . . . . . . . . 30

	II.3.1 Li-ion batteries in electric vehicles . . . . . . . . . . . . . . . . . . . . . . . 30
	II.3.2 Battery thermal management systems (BTMS) . . . . . . . . . . . . . . . . 31
	II.3.3 Physics of immersion cooling . . . . . . . . . . . . . . . . . . . . . . . . . 40

II.4 Chapter conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 II.1 Lithium-ion batteries II.1.1 Generalities and definitions

	A battery is a device capable of converting chemical energy, stored in its active materials, into
	electric energy thanks to electrons transfer generated by oxidation-reduction reactions [89].
	A Lithium-ion battery is a battery producing electricity by a transfer of electrons e -induced
	by chemicals reactions occurring at positive and negative electrodes and involving Lithium
	under the ionic form: the ions Li + , evolving in an electrolyte solution. They are pretty
	light energy storage devices, with a voltage around 4 [V ] and specific energy (expressed in
	watts-hours by kilos) around 100 [W • h • kg -1 ] [124]. They are rechargeable devices, as the
	oxidation-reduction reactions occurring at each electrode are reversible. The electrons and
	Lithium ions transfer is possible in both ways between the electrodes.
	As stated in the introduction, the automotive industry has shown a growing interest in
	these systems for a few years. To widespread the use of Electric Vehicles (EV), the Li-ion
	batteries are required to produce high power during a long time of operation. The electric
	automotive industry seeks to develop vehicles capable of reaching high speeds while widen-
	ing the autonomy range of the vehicles. The capabilities of storage systems in this perspective
	can be compared in Ragone diagrams, like the one in

Table II

 II

					the
	authors sum up the commonly used types of batteries regarding this characteristic. The types
	of batteries are listed with their positive electrode chemistry in Tab. II.1.
	Type of battery	Abbrev. Chemical formula	Voltage	Capacity
				[V]	[W • h • kg -1 ]
	Li. cobalt oxide	LCO	LiCoO 2	3.7-3.9	150
	Li. nickel oxide	LNO	LiN iO 2	3.6	150
	Li. Ni. Co. aluminium oxide NCA	Li(N i 0 , 85Co 0 , 1Al 0 , 05)O 2	3.65	130
	Li. Ni. manganese Co. oxide NMC	Li(N i 0 , 33M n 0 , 33Co o 33)O 2 3.8-4.0	170
	Li. manganese spinel	LMO	LiM n 2 O 4	4.0	120
	Li. iron phosphate	LFP	LiF eP O 4	3.3 / 2.5 -3.6 [19] 130

.1: Types of batteries by positive electrode materials. Extracted from

[START_REF] Wang | A critical review of thermal management models and solutions of lithium-ion batteries for the development of pure electric vehicles[END_REF]

.

  4 gathers the advantages and drawbacks of PCM, indirect liquid, air immersion and liquid immersion

	cooling systems.
	BTMS	Advantages	Drawbacks
	PCM	zero external energy consumption good safety capabilities	low thermal conductivity ensure good ratio thermal conductivities of PCM vs cell
		no secondary cooling loop	low heat transfer
	Air immersion cooling	low cost and design no leak potential	bad temperature homogeneity risk of thermal runaway
		suitable for all type of cells	high space required
			low charging rate
	Indirect liquid cooling	better heat transfer low volume	electric conductive fluid risk of thermal runaway
			bad temperature uniformity
	Liquid immersion cooling	best heat transfer uniform temperature ultra fast charging limits thermal runaway	higher density fluid design complexity require heat sink

Table II.4: Advantages and drawbacks of BTMS. Adapted from

  This first solver, called ICExo, was initially developed by Exoes. The idea behind this model was to get a numerical tool able to give comprehensive insights into the thermal behavior of Li-ion batteries in immersion cooling systems. Also, the evolution of electrical features is of interest to answer Exoes' customer's needs, as they ask for battery packs able to respect some constraints corresponding to specific usage of the vehicles.The solver is seen as a low fidelity (LF) solver as it is based on significant simplifications in solving the conjugate heat transfer between the batteries and the fluid. The electrical parameters of the battery are computed using an Equivalent Circuit Model (ECM, see II.2.1.2). The originality of the present method is that the electrical parameters of the ECM are dependent on pressure and temperature, set as environmental conditions and computed by the model during the operation. The link between these quantities is made thanks to the interpolation of tables provided by battery manufacturers.

III.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 III.2 Thermo-electrical solver: ICExo . . . . . . . . . . . . . . . . . . . . . . . 50 III.2.1 Thermal equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 III.2.2 Electrical equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 III.3 Validation of the model with an experimental test case . . . . . . . . . 58 III.3.1 A review of experimental setups for immersion cooling . . . . . . . . . . . 58 III.3.2 Experimental test case: lab-scale immersion cooling . . . . . . . . . . . . 59 III.3.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 III.3.4 Discussion on the numerical results . . . . . . . . . . . . . . . . . . . . . . 64 III.4 Chapter conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 III.1 Introduction

  is the state of health of the battery, expressed with a percentage. The capacity C is given at the beginning of the simulation by a 1D interpolation in tables linking the temperature and the capacity as follows:C = Interp 1 T ini ; [Temp, C Ah ] .(III.12)Similarly to the 2D interpolation function,Interp1 is based on the 'interp1' function from Matlab. Interp 1 performs a 1D linear interpolation between data in the neighborhood of the query point, by comparison with Interp2, which takes two inputs x, y. The table C Ah is a vector containing values of capacity from 2.88 to 3 [A • h] multiplied by the SOH parameter. The query point is the temperature at t = 0, knowing T ini . The data in the table C Ah are provided by the cell constructor, which one aims to simulate (Murata™ in the case of the experiment presented in the next section).

To compute the voltage of the battery, the Li-ion battery is approximated with an adapted Thevenin model

[58]

. It is represented in Fig.

III

.4 with all the notations of the involved electrical parameters. A parallel RC circuit, with resistance R 1 and capacity C 1 , is connected

  The table R 1 contains values of resistance from 0.5 to 10 [mΩ] with a maximum discrepancy of4 [mΩ] between two values.The query points are the initial temperature T ini and the state of charge SOC(t). Note that the capacity C 1 is the internal capacitance of the battery. It is interpolated from tables linking the temperature, the state of charge and the capacitance, computed as followsC 1 (t) = Interp 2 T ini ,SOC(t); [Temp, SOC, C 1 ] (III.19) The table C 1 contains values of capacitance with a maximum discrepancy of 2.7 • 10 +3 Farads between two values.

	SOC(t); [Temp, SOC, R 1 ]	(III.18)

  Table III.2: Nominal experimental conditions for DS and RC cycles.

  [START_REF] Barbillon | Méthodes d'interpolation à noyaux pour l'approximation de fonctions type boîte noire coûteuses[END_REF]

	Turbulence modeling For some cases, we add a turbulence model to the Navier-Stokes
	equations. The lastly described momentum conservation equation is modified. The Unsteady
	Reynolds Averaged Navier-Stokes (URANS) equations are solved. Under the RANS frame-
	work, the velocity and pressure fields are seen as the sum of a statistically averaged term and
	fluctuations.
	u(x, t) = u(x, t) + u (x, t)	(IV.11)

p(x, t) = p(x, t) + p (x, t) (IV.12)

  rate is imposed as a Dirichlet condition on the air velocity field, as specified in table Tab. IV.3. Then the air exits freely the domain by the boundary ∂Ω outlet . The bottom and top wall boundaries present a no-slip condition. As this computational domain represents only a fraction of a bigger pack, symmetry conditions are applied in the boundaries ∂Ω sides to simulate the interaction with the cells and flow on the sides of this domain.

		Ω f		Ω s
	Ω cylinders coupled condition	Ω cylinders	coupled condition
	∂Ω top	wall -adiabatic	Ω bottom-cylinders adiabatic
	∂Ω bottom wall -adiabatic
	∂Ω sides	symmetry
	∂Ω inlet	u = (0.1; 0; 0) m.s -1
	∂Ω outlet	p = 0
		Table IV.3: Boundary conditions

Table VI .

 VI 1: Nominal experimental conditions for DS and RC cycles.

  Table VI.2: Uncertain parameters and bounds of their respective range of variation.

	Parameter	Symbol Lower bound X min i	Upper bound X max i	Units
	X Mass flow rate	ṁ	2.64 • 10 -3	3.96 • 10 -3	kg • s -1
	X Heat transfer coefficient	h f	200	600	W.m -2 .K -1
	X Internal resistance	R 0	1 • 10 -3	18 • 10 -3	Ω
	X State Of Health	SOH	98	102	%
	X Radial thermal conductivity λ r	0.77	1.55	W.(mK) -1
	X Axial thermal conductivity	λ z	10.37	45.83	W.(mK) -1
	X Specific heat capacity	C p	800	1200	J.(kg.K) -1

  quantities of interest. The total Sobol' indices and first order Sobol' indices for each quantity of interest are shown in Fig. VI.8 and in Fig. VI.9 respectively and detailed in Tab. VI.3 and

	ṁ	h f	R 0	SOH	λ r	λ z	C
	Tab. VI.4.
	ṁ	h f	R 0	SOH	λ r	λ z	C p
	Y 1 0.0000 0.0000 0.0098 0.9879 0.0000 0.0000 0.0000
	Y 2 0.0000 0.0000 0.1360 0.8581 0.0000 0.0000 0.0000
	Y 3 0.1955 0.1668 0.5937 0.0148 0.0029 0.0220 0.0002

Table VI.3: Total Sobol' indices for the three quantities of interest.

The total and first-order Sobol' indices underly the same conclusions about the model's behavior. There is no significant difference between the results from the total and first-order Figure VI.8: Total Sobol' indices for the three quantities of interest. p Y 1 0.0021 0.0021 0.0114 0.9901 0.0021 0.0021 0.0021 Y 2 0.0046 0.0046 0.1383 0.8636 0.0046 0.0046 0.0046 Y 3 0.2053 0.1748 0.5978 0.0236 0.0123 0.0313 0.0096 Table VI.4: First order Sobol' indices for the three quantities of interest.

  Table VII.1: Predicted variance using prior and posterior distributions for some temperatures of interest at given times.

	Output temperature at:	t = 178	t = 377	t = 577	t = 778	t = 1158	t = 1517
	Prior predicted variance	0.001745 0.007332 0.016109 0.027687 0.057450 0.094950
	Posterior predicted variance	0.000077 0.000263 0.000470 0.000667 0.001072 0.001789
	Relative variance reduction [%] 95.6097	96.4173	97.0825	97.5918	98.1342	98.1154

Algorithm 1 Numerical model operation t = 0 Init. thermal and flow parameters:

Compute T (x s , t) through Eq. III.

V.3 Uncertainty forward propagation by sampling

As stated just above, the goal of uncertainty forward propagation is to compute the statistics on the output random variable, and obtain its distribution. The random variable is an outcome of the model M and these statistics could be written as:

Sampling techniques have been widely used to compute these integrals. The advantages of sampling methods are their versatility. However, they present a quite slow convergence rate. More precisely, in the framework of UQ with computer models, the sampling is made from the input random variables.

The step 1, uncertainty characterization, consists in giving a PDF π X i (defining the probability distributions) for each random variable X i of the random vector X = (X 1 , ..., X m ). Then computation of the above integrals is made as follows:

1. Sampling each random variable X i . That is, for each variable X i , draw a number N of samples according to its PDF π X i . This gives a set of values (X (k) i ) k=1,...,N defining the sample set. The sample set of the random vector is similarly written (X (k) ) k=1,...,N .

2. For each sampled input X (k) , evaluate the numerical model. This leads to a sampled set

V.6 Bayesian inverse problem

Let us consider again the numerical model M mapping a m-dimensional input x in a pdimensional output y:

V.6.1 Model calibration

In the framework of model calibration, as seen in the introduction of the chapter, the validation of the numerical model M goes through the comparison of its outputs with respect to some experimental data in a vector y exp . In this work, experimental data is constituted of single observations for each of the p quantities of interest, so y exp ∈ R p . From this perspective, the relation between the model and the observation can be seen as:

ε is called the discrepancy. It accounts for the deviation between the observed measurements and the predictions of the numerical model [START_REF] Wagner | UQLab user manual -Bayesian inversion for model calibration and validation[END_REF], as every model is a simplification of the real physics.

Input uncertainties are represented by the random vector X = (X 1 , ..., X m ) as done earlier in this chapter. The PDF associated with this random vector is denoted π(x). Also, for most real-life applications, the discrepancy ε is not known. A random vector denoted E is used to represent the uncertainties of this additional variable, which is now considered as input of the problem in the sequel. The PDF associated to this random vector is defined by some parameters in a vector x ε and denoted π(x ε ). The covariance matrix associated to this random vector is written K(x ε ).

The PDFs π(x) and π(x ε ) represent the information known about the unknowns, before having observed the data y exp . They are called the prior distributions. Also, in the Bayesian paradigm, the experimental data are realizations of the random variable Y exp representing the data in y exp .

We denote x I the vector containing the realizations of the model inputs and the discrepancy variable, and X I the associated random vector:

The goal of the Bayesian calibration is to find the probability of the input variables X I , conditioned to the information provided by Y exp . In other words, it is to find the PDF associated with the random variable:

CALIBRATION OF THERMO -ELECTRICAL PARAMETERS UNDER UNCERTAINTIES

Chapter abstract

In this chapter, we apply the LF numerical model on the experimental test case from Exoes, presented in the chapter III. Uncertainties coming from the input multi-physics parameters of the model are considered. First, we assess the impact of the uncertainties on the temperature evolution by considering the prior input distributions. Then, we directly perform a Bayesian calibration of these parameters using the experimental measurements. The informative distributions outcoming of this process are used to validate the model on the two experimental conditions available. The uncertainty in the model's temperatures prediction is significantly reduced from this process. Finally, for a pure design purpose, the learned distributions of inputs and the numerical model are used to simulate the system under conditions representing a realistic racing car operation. A Sobol indices based sensitivity analysis is performed to get further analysis elements on the behavior of the BTMS. The results presented in this chapter have been submitted for publication in a journal.

Outline

VI.2.5 Global sensitivity analysis

Specifically, in the framework of the sensitivity analysis study, we are not interested anymore in vector outputs of the black box model but scalar outputs. Then, in this part, the black-box model is a function mapping the seven inputs in X into a scalar quantity. The according model is denoted M Y . Considering those notations, the ANOVA decomposition of the model M Y can be written as follows:

second order

seventh order .

(VI.9)

Let's write here the definition of the Sobol indices seen in chapter V, with the notations of this section. The first order Sobol index S 1 i for an input variable X i is defined by:

(VI.10)

The total Sobol index S T i for the variable X i is computed as follows:

where X\i represents the set of all input variables and their combinations except X i .

VI.3 Results

VI.3.1 Posterior distributions calibrated from the Datasheet cycle case.

The Bayesian inverse problem is solved using the Datasheet cycle case (Tab. VI.1) which the experimental results are described in III.3.2. In this section, the simulation is performed only on the first 250 seconds of the cycle during the discharging sequence. The experimental data T exp selected to learn the distributions is constituted of six measurements of temperatures at different times t exp = [START_REF] Goodman | Ensemble samplers with affine invariance[END_REF][START_REF]MCMC algorithms[END_REF][START_REF] Stroe | A comprehensive study on the degradation of lithium-ion batteries during calendar ageing: The internal resistance increase[END_REF]159,199,239] seconds. For computational reasons related to the performances of the MCMC algorithm, the choice has been made to compute the posterior distributions with only six measurements of the discharging sequence. In The first consideration from these results is that the supports of the posterior probability density functions are significantly more narrow than the prior uncertainties arbitrarily defined in

INTERNAL RESISTANCE MODEL CALIBRATION UNDER

UNCERTAINTIES

Chapter abstract

This chapter is devoted to the application of UQ techniques to the test case presented in section IV.4 with the CFD solver. We propose an innovative way to deal with the uncertainties related to internal resistance. A model for the internal resistance dependent on temperature is constructed, allowing the treatment of the uncertainties coming from this parameter. Prior uncertainties on the inputs are first propagated thanks to a surrogate model. The prior distributions lead to a considerable variation in the temperature prediction. Then, Bayesian calibration of the internal resistance model parameters using experimental temperature measurements is performed. The uncertainty in the inputs is significantly reduced in the outcoming posterior distributions. These posterior distributions are propagated, and we obtain a much more narrow confidence interval in the temperature predicted by the numerical model. Finally, an exploratory action is conducted to construct an internal resistance model dependent on both the state of charge and the temperature. This model is implemented in the CFD solver, and we can assess the impact of low SOC on temperature evolution. The results presented in this chapter have been submitted for publication in a journal. We observe in Section IV.4 that the 3D configuration permits to represent the effects due to natural convection, contrary to the 2D. On the contrary, we expect that the temperature prediction will be more influenced by the wide range of values taken by the resistance model, due to uncertainties, than the slope changes caused by the natural convection represented in 3D. Furthermore, as stated earlier in this manuscript, the uncertainty forward propagation approach requires numerous runs of the CFD solver. Running the simulations in 2D represents a significant reduction in computational costs than in 3D. The construction of the Design of Experiment with the 3D simulations was then considered too costly. Finally, the

Outline

VII.2.2 Construction of the internal resistance model

Let us now focus on the construction of the internal resistance model, which is input of the source term for the heat equation in the solid domain. As described earlier, we recall that the source term represents the ohmic heating of the batteries through the internal resistance and input electric currents evaluations:

So here, we want to construct a mathematical function representing the internal resistance R. When dealing with experiments not reaching very low values of state of charge for the batteries, the sole temperature dependence stands as a reasonable assumption. Several references studying properties of Li-ion batteries assume that the internal resistance is dependent on the temperature, in a cubic polynomial form [START_REF] Park | Dynamic Thermal Model of Li-Ion Battery for Predictive Behavior in Hybrid and Fuel Cell Vehicles[END_REF][START_REF] Jilte | Numerical Investigation on cooling performance of Li-ion battery thermal management system at high galvanostatic discharge[END_REF][START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF]. The objective here is then to build a direct relationship between the resistance of the cell and its temperature, following some physical constraints, i.e. T → R(T ). Specifically, the built model R(T ) is imposed to be strictly decreasing within the temperature range of the experiment following several works [START_REF] Park | Dynamic Thermal Model of Li-Ion Battery for Predictive Behavior in Hybrid and Fuel Cell Vehicles[END_REF][START_REF] Karimi | Thermal management of lithium-ion batteries for electric vehicles[END_REF][START_REF] Li | Effects of temperature on dynamic characteristics of li-ion batteries in electric vehicle applications[END_REF][START_REF] He | Combined experimental and numerical study of thermal management of battery module consisting of multiple Li-ion cells[END_REF]. These works constructed their model from direct internal resistance measurements at different temperatures.

A Bezier parametrization is chosen to model the resistance, requiring the definition of some control points. Note that the expected monotonic behavior of the resistance can be naturally imposed with a suitable choice of the control points. Within the range of rele-

VII.2.3 Uncertainty representation through the internal resistance model

Uncertainties of the internal resistance can be represented by modelling the distribution of the control points variability using the Bezier parametrization. Using the probabilistic framework defined in chapter V, a random variable R i is introduced, giving the resistance value at the temperature abscissa T i , i ∈ [0, 3] of each control point.

The random variables R i map the event ω from a probability space into a value of resistance in [mOhm]: R i : ω → R i (ω). Then we construct the random vector composed by the four resistance values:

The uncertainty coming from the internal resistance parameter is then represented through are quite variable. Then the source term evaluations present significant differences from one sample run to another. These differences are traduced in the temperature prediction by significant slopes changes between all the samples. As time increases, the range of predicted temperatures gets bigger.

So the current results could be not satisfying for a robust design, as the amount of information provided by the prior distributions increases the uncertainty significantly on the temperature with the simulated time. The following section performs the calibration process to maintain this uncertainty in a smaller range over time.

VII.4 Calibration of the resistance model using experimental data

In this section, the approach considered consists of using the experimental data directly to reduce this variability in the temperature prediction. We illustrate two methods to incorporate the information brought by the experimental measurements to correct the temperature prediction. The difference between a deterministic and Bayesian calibration is shown, from the implementation to the results and resulting analysis that can be done accordingly.

VII.5 Effect of SOC on the Li-ion cells temperature evolution

So far, only the temperature has been considered, which is the input variable of the internal resistance. This assumption yields convincing results in predicting temperature for experiments where a very low state of charge (SOC) is not reached. Although, the state of charge might still significantly affect the resistance value when performing full discharge cycles. This impact on the resistance value is expected to change the shape of the temporal temperature response from the heated Li-ion batteries. Many studies have proposed methods or experimental investigations to represent the effect of SOC in the model for the internal resistance [START_REF] Li | Effects of temperature on dynamic characteristics of li-ion batteries in electric vehicle applications[END_REF][START_REF] Ansean | Electric Vehicle Li-Ion Battery Evaluation Based on Internal Resistance Analysis[END_REF][START_REF] Wu | Improved State of Charge Estimation for High Power Lithium Ion Batteries Considering Current Dependence of Internal Resistance[END_REF][START_REF] Yang | Lithium-Ion Battery Internal Resistance Model Based on the Porous Electrode Theory[END_REF][START_REF] Stroe | A comprehensive study on the degradation of lithium-ion batteries during calendar ageing: The internal resistance increase[END_REF][START_REF] Chen | Estimation the internal resistance of lithium-ion-battery using a multi-factor dynamic internal resistance model with an error compensation strategy[END_REF]. The techniques used in these references mainly measure the values of resistance at different SOC and temperature points. However, examining these references, it is hard to claim a general behavior of the internal resistance depending on the SOC dependency. Some studies show decreasing trends, as others illustrate a strong dependency but more complex relationships.

This section proposes a method to consider the effect of SOC and build an internal resistance model accordingly, also based on experimental measurements. The present method makes easy the implementation of the constructed model in the CFD solver and allows to predict the temperature response of the heated batteries.

VII.5.1 Construction of an internal resistance model from experimental data

Unlike the previous methodology, where we tried to reconstruct the shape of the internal resistance model indirectly from temperature measurements, here we start from internal resistance measurements. We try to construct an internal resistance model accounting for the effect of SOC and temperature based on direct measurements on internal resistance depending on those two parameters. In this part, we consider a pure deterministic framework, where the model's coefficients are not resulting anymore from a stochastic process.

Here the Design of Experiment is the set containing these experimental points. The Kriging is not used to replace a computational model but as a pure interpolation tool to construct a model based on given data.

We obtain a Kriging predictor R which depends on the inputs z = (SOC, T ). The parametrization of this function relies on the definition of correlation functions, describing the correlated effect of two input points z and z on the output R. Then, the construction of the Kriging model consists of optimizing the coefficients (the hyperparameters) by maximizing the likelihood of observing the construction points with the model. The Kriging predictor is a random variable indexed by z following a Gaussian distribution of mean µ R(z) and vari-

Once the Kriging model is constructed, it is evaluated on a large set of test points X test .

This set contains N = 1 • 10 4 points, sampled following a uniform distribution over the range of variations considered for each input variable. In In terms of computational and implementation efforts, it is not convenient to input the Kriging predictor directly in the CFD model. Indeed, supplementary matrix operations would be necessary to evaluate the Kriging function R(SOC, T ) and compute the source term in the heat equation at each time step. For this reason, the surface generated by the Kriging estimator on the N points from X test is fitted with bi-variate polynomials, using least square regression. The polynomial structure is chosen a priori: a maximal degree 5 is considered for the SOC dependency, and degree 3 for the temperature. The bi-variate polynomial fitting the Kriging predictor surface is denoted P 53 . The quality of the fit is assessed, and the final

CONCLUSION AND PERSPECTIVES

Chapter abstract

This thesis addressed the issue of numerical simulation under uncertainties of immersion cooling systems for Li-ion batteries. This concluding chapter covers at a glance the contributions of this work. Then, some perspectives are given to illustrate some potential future developments.

Outline

Numerical simulation for the immersion cooling of Li-ion batteries

The first objective of the thesis was to develop numerical tools to address the simulation of immersion cooling systems of Li-ion batteries. Two main contributions have allowed addressing this issue.

Fast and low fidelity solver for immersion cooling (chapter III)

First, in close collaboration with Exoes, a fast and accurate numerical model (ICExo) has been developed to simulate the heat transfer for Li-ion batteries immersed in a cooling fluid. The main novelty of this model is its original approach. The modeling strategy for the heat transfer resolution significantly reduced computational costs compared to a complete conjugate heat transfer approach. Here, the heat equation is solved in the solid domain, and the heat transfer between the fluid and the solid is computed through enthalpy balances along a conduct. The solver still presents a good ability to reproduce experimental results in terms of temperature prediction. Also, this model can simulate some electrical parameters evolution (voltage and state of charge).

Practically, the versatility of this numerical solver allowed us to study the thermal and electrical behavior of the simulated battery pack for complex input conditions, i.e. unsteady input electrical currents representing realistic operation conditions of electric vehicles.

This low fidelity solver gives a good overview of the general operation of a battery pack immersed in a flow of cooling fluid in terms of thermal and electrical performances. In that the battery's internal resistance and state of health. The thermal parameters influence the quantity related to temperature: mass flow rate, heat transfer coefficient, thermal conductivities, and specific heat but with less importance. The most critical parameter is the internal resistance, which governs the heat source produced in the batteries.

Bayesian calibration of the multi-physics parameters and uncertainty reduction (chapter VI) In parallel to the uncertainty propagation approach, the original experimental data from Exoes allowed calibrating the low fidelity solver while considering the uncertainties on the input parameters. Prior uncertainties have been assumed for the input parameters, based on literature review and the knowledge brought by Exoes experts. The distributions characterizing these uncertainties were quite large and gave model responses with a large variability (see the previous contribution). The Bayesian calibration process gave more narrow and informative distributions on the inputs values. Also, the forward propagation of these posterior distributions significantly reduced the width of the 95% confidence intervals of temperature prediction. These distributions were also used on different physical conditions, and the model's predictions showed a good comparison with the experimental measurements after the calibration process. Finally, we showed that the information brought by the experimental data allowed us to reduce the uncertainty on the inputs parameters and the temperature predicted by the low fidelity solver. Note that the sensitivity analysis previously stated was performed using these distributions, giving much more informative results on the estimation of the Sobol indices.

Accuracy assessment of an internal resistance model (chapter VII) Finally, using the CFD solver for the conjugate heat transfer problem, we proposed a method to construct an internal resistance model involved in the source term of the heat equation. The methodology allowed to build an internal resistance model dependent on temperature. The Bezier parameterization eased the control of its shape to follow physical constraints. Uncertainties related to the resistance parameter were represented using this model construction. Prior uncertainties were assumed for each model parameter based on the CFD solver's temperature response.

Uncertainties were propagated using sampling methods through a surrogate model. The 95% confidence interval with prior uncertainties embraced the experimental measurements and error envelopes well. However, as the simulation is transient, the uncertainty traduced by the width of the confidence interval showed a significant increase over time.

A Bayesian calibration was performed to reduce the prior uncertainties assumed on the parameters of the resistance model. Based on the experimental measurements of temperature, the posterior distributions were computed. The posterior distributions of the resistance values showed a significant variance decrease compared to the priors. The shape of the distributions also showed a much more peaky shape. Also, the propagation of the posterior distributions gave a 95% confidence interval on the temperature prediction with a quasisteady width over time while embracing the experimental data systematically. Experimental measurements allowed gaining significant information on the behavior of the resistance parameter and the resulting temperature prediction.

VIII.2 Perspectives

We draw here some perspectives for future works.

Calibration of internal resistance model applied on an industrial-scale immersion cooling problem

The process performed for the calibration of the internal resistance model (chapter VII) can be done on a test case closer to a real industrial application. Exoes expressed some needs in the simulation of natural convection for batteries set up in a realistic module of immersion cooling. TrioCFD showed good abilities to reproduce the natural convection for such configurations (section IV.4).

The idea of this contribution would be to gather the experimental data of temperature provided by Exoes from this module, the UQ framework of Bayesian calibration for the internal resistance model developed during this thesis, and the CFD solver applied on this geometry. Then we would be able to calibrate the internal resistance model parameters using experimental data and the simulations, considering the uncertainties on the resistance input parameters. Finally, after the calibration, we can generate a robust temperature prediction.

Then, one can assess if the CFD predictions, including the reduced uncertainty, match the industrial constraints the BTMS has to fulfill.

Calibration of internal resistance model including the SOC effect

In the section VII.5

(chapter VII), we proposed a method to include the dependence of the state of charge in the internal resistance model and assess the effect of this parameter on the temperature evolution of the heated batteries. This enhanced resistance model in the source term was set up as a bi-variate polynomial to take account of the dependency from temperature and state of charge. This enhanced internal model could also be calibrated, solving the Bayesian inverse problem, similarly to Chapter VII. Obtaining meaningful results from calibration requires a full discharging sequence of the batteries during the experimental measurements. In this way, we could observe the effect of the low state of charge on the resistance values and then on the resulting temperature response. We recall that in the present work, the experimental temperature measurements were available for a discharging sequence going down only to a state of charge of 30%. Yet, experimental data of resistance supporting the construction of this enhanced model showed that the effect of state of charge gets significant under 30%. We can expect that the Bayesian calibration process would also lead to narrow distributions for all the resistance control points parametrizing the bivariate polynomial. This anisotropic operator has been implemented in TrioCFD by collaborators at CEA. potential contribution would be to merge the information provided by the two models. The low fidelity solver presents a small computational cost compared to the high fidelity solver.

Anisotropic thermal conductivity of

Then, merging the predictions of both models could be interesting in terms of computational resources while keeping a high level of fidelity overall.

In practice, this would be implemented using the same approach as the work presented in [START_REF] Le | Multi-fidelity Gaussian process regression for computer experiments[END_REF]. Firstly, the idea is to set up both models as a black box, taking the same inputs and producing the same nature of outputs. Here, the approach performed with the CFD solver in chapter VII would be reproduced with the low fidelity solver. The inputs are the four resistance values parametrizing the internal resistance model, and the outputs are a vector of temperatures evaluated at different times of the transient simulation.

Then, we build two designs of experiments using both solvers. The Design of Experiments built with the high fidelity solver will naturally contain fewer points than the Design of Experiments of the low fidelity solver to keep a good computational cost accuracy trade-off.

Then a Kriging surrogate model can be built, merging the information learned from those two Design of Experiments. Indeed, as presented in [START_REF] Lataniotis | The Gaussian process modelling module in UQLab[END_REF], a so-called hierarchical Kriging can be setup. First, we build a Kriging model from the evaluations of the low fidelity solver. Then, another Kriging model is constructed using the evaluations of the high fidelity solver, but its

Résumé

Pour promouvoir une utilisation plus large des véhicules électriques, les batteries Lithium-ion (Li-ion) se doivent de supporter des courants électriques importants, générant ainsi de fortes contraintes thermiques qui dégradent leurs performances et leur durée de vie. La gestion thermique des packs de batteries est donc un élément crucial pour répondre à ces nouvelles contraintes industrielles. La technologie de refroidissement par immersion est une solution prometteuse en termes de performances thermiques. Ces systèmes sont gouvernés par des phénomènes multi-physiques, allant de la chimie interne des batteries jusqu'au transfert thermique à l'échelle du pack de batteries alimentant le moteur électrique. Cette thèse a pour objectif de développer des modèles numériques pour le refroidissement par immersion des batteries Li-ion, tout en considérant les incertitudes provenant des paramètres physiques en jeu.

Cette problématique est abordée en proposant deux modèles de fidélités croissantes. D'abord, un modèle dit basse fidélité est développé, incluant la modélisation des phénomènes thermiques et électriques du problème. Des méthodes de quantification d'incertitudes (calibration Bayésienne et analyse de sensibilité) couplées avec des données expérimentales originales offrent ainsi des éléments de compréhension et d'analyse sur le comportement global du système. Ensuite, une approche plus spécifique est présentée à l'aide d'un code de calcul CFD haute fidélité. Le calcul du transfert thermique conjugué sous régime transitoire en deux dimensions d'un pack de batteries immergées est ainsi réalisé. Cet outil est utilisé pour évaluer la précision d'un modèle construit a priori, représentant la résistance interne des batteries Li-ion. Les incertitudes provenant de la résistance sont prises en compte grâce à la paramétrisation de ce modèle et calibrées en utilisant un cas test expérimental de la littérature. Enfin, pour obtenir une meilleure compréhension de la physique de ces problèmes de refroidissement par immersion, la fidélité de l'outil CFD est augmentée en considérant des calculs de transfert thermiques en 3D, ainsi qu'un modèle de résistance interne amélioré.

Mots clés Batteries Lithium-ion; Refroidissement par immersion; Simulation numérique; Quantification d'incertitudes; Transfert thermique conjugué; Calibration Bayésienne; Méta-modèles.

NUMERICAL SIMULATION AND UNCERTAINTY QUANTIFICATION FOR IMMERSION COOLING OF LITHIUM-ION BATTERIES

Abstract

To encourage a wider use of electric vehicles, Lithium-ion (Li-ion) batteries are required to handle high electric currents, generating great heat loads which deteriorate their performances and lifespan. The thermal management of the battery packs is a key element to fulfill these industrial demands. Immersion cooling technology stands as a promising solution in terms of heat transfer performances. Multi-physics processes govern those systems, from the internal chemistry of Li-ion cells to the heat transfer at the battery pack scale powering the electric engine. This thesis aims to develop numerical models of immersion cooling systems for Li-ion batteries considering the uncertainties coming from the physical parameters.

This issue is addressed by proposing two models of increasing fidelity. Firstly, a low fidelity model is developed, including the thermal and electrical phenomena of the immersion cooling problem. Uncertainty quantification methods (Bayesian calibration and sensitivity analysis) coupled with original experimental data provide a deeper knowledge on the overall behavior of the system. Secondly, a more specific approach is performed using a high fidelity Computational Fluid Dynamics (CFD) model solving the transient conjugate heat transfer in an immersed battery pack in two dimensions. This CFD tool is used to assess the accuracy of a constructed model for the internal resistance of Li-ion batteries. Uncertainties coming from the internal resistance are taken into account thanks to the parameterization of this model and calibrated using an experimental test case from literature. Furthermore, for a better understanding of the immersion cooling physics, the fidelity of the CFD model is increased by considering 3D simulations and an enhanced internal resistance model.

Keywords Lithium-ion batteries; Immersion cooling; Numerical simulation; Conjugate heat transfer; Uncertainty quantification; Bayesian calibration; Surrogate modeling.

Inria, Institut de Mathématiques de Bordeaux

UMR 5251 Université de Bordeaux, 33400 Talence, France.