Thermal management study of an electrolytic stack using computational fluid dynamics

Master Thesis Master’s degree in Industrial Engineering Thermal management study of an electrolytic stack using computational fluid dynamics REPORT January 25, 2024 Autor: Antonio Roig Andreu UPC Tutor: Elisabet Mas de les Valls Delivered: January 2024 Escola Tècnica Superior d’Enginyeria Industrial de Barcelona Abstract This master’s thesis has been carried out in relation to the studies of the master’s degree in industrial engineering. The project comes from an interest in understanding how electrolyzers work to produce hydrogen and to develop a model that analyzes, by means of computational fluid dynamics (CFD), the heat management of these devices. Therefore, it will first be necessary to study the physical and chemical fundamentals of electrolysis and and the different types of water electrolyzers for hydrogen production that exist nowadays. In addition, since the study focuses on the use of CFD for thermal analysis, it will be explained what the physical fundamentals of computational fluid dynamics are, and how the computer programs that allow these fundamentals to be applied work. In particular, the project focuses on a SOEC type electrolyzer. In order to model these devices, it has been necessary to understand which layers of materials make up these devices, as well as the function that each of them performs in the electrolysis of water. Some of the thermodynamic and electrical properties of the different layers have been obtained as a function of temperature. Subsequently, the different SOEC configurations that currently exist have been analysed, and the geometry to be modelled has been defined. The physical-theoretical equations that allow the behaviour of these devices to be simulated have also been analyzed. These equations have been explained in detail. They have then been simplified in order to be able to apply them to the model that has been necessary to carry out the simulation of this project and, in this way, to obtain a first approximation of the thermal analysis of this type of electrolyzers. In order to check whether the modelling idea proposed is correct, a base model has first been developed to verify the correct use of the tools provided by the OpenFOAM software. In this base model, three solid regions emulating the SOEC regions have been defined. The meshing and how the simulation has been prepared has been explained, and then temperature and voltage results have been obtained. The same concepts applied in this base model were then used for the development of the SOEC model. In this model, five regions have been defined, with the corresponding transport equations: the two bipolar plates, the membrane-electrolyte assembly, the supplied water channels and the produced hydrogen channel. Finally, points for improvement of the model have been identified and future implementations have been proposed, which could lead to more accurate results showing the real thermal behaviour of these devices. page. 2 Report Resumen La presente tésis de máster se ha llevado a cabo en relación con los estudios del máster de ingeniería industrial. El proyecto parte del interés por comprender cómo funcionan los electrolizadores para la producción de hidrógeno y por la posibilidad de desarrollar un modelo que analice, mediante fluido-dinámica computacional (CFD), el comportamiento térmico de estos dispositivos. Por lo tanto, primero será necesario estudiar los fundamentos físico-químicos y los distintos tipos de electrolizadores de agua para la producción de hidrógeno. Además, puesto que el estudio se centra en la utilización de CFD para el análisis térmico, se explicará cuales son los fundamentos físicos de la fluido-dinámica computacional, y cómo trabajan los programas informáticos que permiten aplicar dichos fundamentos. Concretamente, el proyecto se centra en un electrolizador de tipo SOEC. Para poder modelar dichos dispositivos ha sido necesario entender qué capas de materiales componen estos dispositivos, así como la función que ejercen en la electrólisis del agua cada uno de estos. Se han obtenido algunas de las propiedades termodinámicas y eléctricas de las distintas capas en función de la temperatura. Posteriormente, se han analizado las diferentes configuraciones de SOEC que existen actualmente, y se ha definido la geometría a modelar. También se han analizado las ecuaciones físico-teóricas que permiten simular el comportamiento de estos dispositivos. Estas ecuaciones se han explicado de manera detallada. A continuación, se han simplificado con el objetivo de poder aplicarlas en el modelo que ha sido necesario realizar para poder llevar a cabo la simulación de este proyecto, y de esta manera, conseguir una primera aproximación del análisis térmico de este tipo de electrolizadores. Para comprobar si la idea de modelado planteada es correcta, se ha desarrollado primero un modelo base que ha permitido verificar la correcta utilización de las herramientas que el software OpenFOAM facilita. En este modelo base se han definido tres regiones sólidas que emulan las regiones del SOEC. Se ha explicado el mallado y cómo la simulación ha sido preparada, y posteriormente se han obtenido unos resultados de temperatura y voltaje. Los mismos conceptos aplicados en este modelo base han sido utilizados después para el desarrollo del modelo del SOEC. En dicho modelo se han definido cinco regiones distintas con las correspondientes ecuaciones de transporte: las dos placas bipolares, el ensamble membrana-electrolito, los canales de agua suministrada y el canal de hidrógeno producido. Finalmente, se han identificado puntos de mejora del modelo y se han propuesto implementaciones para llevar a cabo en un futuro, las cuales podrían permitir obtener resultados más precisos que muestren el comportamiento térmico de estos dispositivos. Thermal management study of an electrolytic stack using CFD page. 3 Resum Aquesta tesi de màster s’ha dut a terme en relació amb els estudis del màster d’enginyeria industrial. El projecte part de l’interès per comprendre com funcionen els electrolitzadors per a la producció d’hidrogen i per la possibilitat de desenvolupar un model que analitzi, mitjançant fluid-dinàmica computacional (CFD), el comportament tèrmic d’aquests dispositius. Per tant, primer serà necessari estudiar els fonaments fisico-químics de l’electròlisi i els diferents tipus de electrolitzadors d’aigua per a la producció d’hidrogen. A més, ja que l’estudi se centra en la utilització de CFD per a l’anàlisi tèrmica, s’explicarà quals són els fonaments físics de la fluid-dinàmica computacional, i com treballen els programes informàtics que permeten aplicar aquests fonaments. Concretament, el projecte es focalitza en un electrolitzador de tipus SOEC. Per a poder modelar aquests dispositius ha estat necessari entendre quines capes de materials composen aquests dispositius, així com la funció que exerceixen en l’electròlisi de l’aigua cadascun d’aquests. S’han obtingut algunes de les propietats termodinàmiques i elèctriques de les diferents capes en funció de la temperatura. Posteriorment, s’han analitzat les diferents configuracions de SOEC que existeixen actualment, i s’ha definit la geometria a modelar. També s’han analitzat les equacions físic-teòriques que permeten simular el comportament d’aquests dispositius. Aquestes equacions s’han explicat de manera detallada. A continuació, s’han simplificat amb l’objectiu de poder aplicar-les en el model que ha estat necessari realitzar per a poder dur a terme la simulació d’aquest projecte, i d’aquesta manera, aconseguir una primera aproximació de l’anàlisi tèrmic d’aquesta mena d’electrolitzadors. Per a comprovar si la idea de modelatge plantejada és correcta, s’ha desenvolupat primer un model basic que ha permès verificar la correcta utilització de les eines que el programa OpenFOAM facilita. En aquest model basic s’han definit tres regions sòlides que emulen les regions del SOEC. S’ha explicat l’emmallat i com la simulació ha estat preparada, i posteriorment s’han obtingut uns resultats de temperatura i voltatge. Els mateixos conceptes aplicats en aquest model basic han estat utilitzats després per al desenvolupament del model del SOEC. En aquest model s’han definit cinc regions amb les corresponents equacions de transport: les dues plaques bipolars, l’acobli membrana-electròlit, els canals d’aigua subministrada i el canal d’hidrogen produït. Finalment, s’han identificat punts de millora del model i s’han proposat implementacions per a dur a terme en un futur, les quals podrien permetre obtenir resultats més precisos que mostrin el comportament tèrmic real d’aquests dispositius. page. 4 Report Contents 1 Preface 10 1.1 Origin of the project .................................... 10 1.2 Motivation ......................................... 10 2 Introduction 12 2.1 Goals/Objectives ..................................... 12 2.2 Scope of the project .................................... 12 2.3 State of the art ....................................... 13 3 Background 15 3.1 Water electrolysis to produce hydrogen ........................ 15 3.1.1 Alkaline electrolyzer ............................... 15 3.1.2 Proton Exchange Membrane electrolyzer ................... 17 3.1.3 Solid Oxide Electrolyzer ............................. 19 3.2 OpenFOAM ........................................ 22 3.2.1 Fundamentals ................................... 22 3.2.2 Case structure ................................... 24 3.2.2.1 Simulation constants (constant folder) ............... 25 3.2.2.2 Boundary conditions (0 folder).................... 26 3.2.2.3 Simulation parameters (system folder) ............... 26 3.2.3 Mesh ........................................ 27 3.2.4 Post-processing .................................. 29 4 SOEC electrolyzer model 30 4.1 Layers and materials ................................... 30 4.1.1 The hydrogen electrode (Cathode) ....................... 31 4.1.2 Electrolyte ..................................... 34 4.1.3 The intermediate layer .............................. 36 4.1.4 The oxygen electrode (Cathode) ........................ 38 4.1.5 The interconnector ................................ 40 4.2 Geometry of the model .................................. 43 4.3 Equations of the model .................................. 45 4.3.1 Equations of the theoretical model without simplifying ........... 45 4.3.2 Simplified theoretical model equations .................... 52 5 Base model 57 5.1 Definition of the BM ................................... 57 5 page. 6 Report 5.2 Meshing the BM ...................................... 58 5.3 Preparation of the BM simulation ............................ 64 5.4 Simulation results for BM ................................ 70 6 SOEC model 73 6.1 Definition of the SOEC model .............................. 73 6.2 Meshing the SOEC model ................................ 74 6.3 Preparation of the SOEC model simulation ...................... 77 6.4 Simulation results for SOEC ............................... 81 6.5 Future steps ........................................ 84 7 Economic assessment 86 7.1 Human resources cost .................................. 86 7.2 Energy cost ......................................... 86 8 Environmental impact 88 9 Social and gender equality impact 89 10 Project planning 90 Conclusions 91 Acknowledgements 92 References 93 List of Figures 1 Schematic representation of the operation of an alkaline electrolysis cell [8]. . . . 16 2 Schematic representation of the configuration and operation of a PEM electrolysis cell [10]. .......................................... 18 3 Schematic representation of the configuration and operation of a SOEC electrolysis cell [12]. ........................................ 20 4 Structure for the preparation of a CFD simulation in OpenFOAM [16]. . . . . . . 27 5 A single block generated using blockMesh esh...................... 28 6 Meshing process using snappyHexMesh......................... 28 7 Scanning electron microscopy (SEM) cross section image of the different layers constituting a cathode-supported cell, from bottom to top: Ni-YSZ, hydrogen, electrode (support + functional), YSZ electrolyte, YDC intermediary layer, and LSCF oxygen electrode [12]. ............................... 30 8 Density in function of temperature for YSZ-Ni. .................... 32 9 Specific heat in function of temperature for YSZ-Ni. ................. 33 10 Thermal conductivity in function of temperature for YSZ-Ni. ............ 33 11 Electrical resistivity in function of temperature for YSZ-Ni. ............. 33 12 Density in function of temperature for YSZ. ...................... 35 13 Specific heat in function of temperature for YSZ. ................... 35 14 Thermal conductivity in function of temperature for YSZ. .............. 35 15 Electrical resistivity in function of temperature for YSZ. ............... 36 16 Density in function of temperature for GDC. ..................... 37 17 Specific heat in function of temperature for GDC. ................... 37 18 Thermal conductivity in function of temperature for GDC. ............. 38 19 Electrical resistivity in function of temperature for GDC. .............. 38 20 Specific heat in function of temperature for Cobalt [22]. ............... 39 21 Thermal conductivity in function of temperature for Cobalt [22]. ......... 40 22 Electrical resistivity in function of temperature for Cobalt [22]. ........... 40 23 Density in function of temperature for Stainless Steel. ................ 41 24 Specific heat capacity in function of temperature for Stainless Steel. ........ 42 25 Thermal conductivity in function of temperature for Stainless Steel. ........ 42 26 Electrical resistivity in function of temperature for Stainless Steel. ......... 42 27 Geometry of the SOEC that will be modeled [5]. ................... 43 28 Detail of the layers that constitude the membrane. .................. 44 29 BM definition for simulation. .............................. 58 30 STL file generated for the BM. .............................. 59 31 blockMesh generated for the BM. ............................ 61 7 32 blockMesh and STL file overlaid for BM. ........................ 61 33 Mesh generated with snappyHexMesh for the BM..................... 63 34 Regions created with splitMeshRegions for the BM..................... 64 35 Results obtained for temperature in BM. ........................ 71 36 Results obtained for voltage in BM. ........................... 72 37 STL file generated for the BM. .............................. 74 38 blockMesh generated for the SOEC model. ....................... 75 39 blockMesh and STL file overlaid for BM. ........................ 75 40 Mesh generated with snappyHexMesh for the SOEC. ................. 76 41 Regions created with splitMeshRegions for the SOEC. ................. 76 42 Temperature results obtained for SOEC. ........................ 82 43 Temperature map of the SOEC cross section. ..................... 83 44 Voltage results obtained for SOEC. ........................... 83 45 Voltage map of the SOEC cross section. ........................ 84 46 Gantt diagram for the project. .............................. 90 List of Tables 1 SOEC dimensions to be modelled ............................ 44 2 BM geometry and operating parameters. ........................ 58 3 BM thermophysical properties for each region. .................... 65 4 Boundary conditions for SOEC. ............................. 74 5 Thermophysical properties for each MEA layer. .................... 77 6 SOEC thermophysical properties for each region. ................... 78 7 Electrical resistivities for each MEA layer. ....................... 79 8 Human resources cost summary ............................ 86 9 Energy cost summary ................................... 86 10 Final costs ......................................... 87 11 Results for CO2emissions and nuclear waste produced ............... 88 8 Thermal management study of an electrolytic stack using CFD page. 15 3 Background 3.1 Water electrolysis to produce hydrogen Hydrogen electrolysis is a chemical process that uses electricity to produce a non-spontaneous chemical reaction. In this process, water is splited into hydrogen and oxygen gases, by applying an electrical current to the anode and cathode side. It is important to note that, depending on the technology used to produce the electrolysis of hydrogen, the reactions that occur at the cathode and anode of the cell vary. However, the general reaction that takes place in any water electrolytic process is shown in equation 1. 2H2O+elecrical energy →2H2+O2(1) The device where this reaction takes place is called electrolyzer. In general, it consists of an anode and a cathode that are separated by an electrolyte, which can be a liquid or a solid material. The basic operation consists on applying an electric current to the electrodes, which causes the water molecules to split into their constituent elements, producing in the cathode electrode hydrogen gas and oxygen gas in the anode electrode. Nowadays, there are three main types of electrolizers which are as follows [6]: •Alkaline Electrolyzers: In alkaline electrolyzers an alkaline electrolyte solution (normally potaisum) is used as the electrolyte. It operates via transport of hydrogen ions (OH−) through the electrolyte. More information about this kind of electrolyzers can be found in section 3.1.1. •Proton Exchange Membrane Electrolyzers: This kind of electrolyzers uses a proton exchange membrane as the electrolyte. In this case, it takes place the transport of protons (H+) from the anode to the cathode through the electrolyte that separates them. More information about this kind of electrolyzers can be found in section 3.1.2. •Solid Oxide Electrolyzers: In this type pf technology, a solid ceramic material is used as the electrolyte, which conducts negatively charged oxygen ions (O2−) from the cathode to the anode. More information about this kind of electrolyzers can be found in section 3.1.3. 3.1.1 Alkaline electrolyzer Alkaline water electrolyzers are the most ancient existing devices to produce the electrolysis of water. It is a really mature technology that produces high-purity hydrogen, and nowadays, is page. 16 Report the most widely applied today to produce it in industrial processes [7]. Like any electrolytic cell, the alkaline electrolysis cell consists of two electrodes (an anode and a cathode) separated by an electrolyte. Direct current is supplied to these electrodes to produce the electrolysis of water. The particularity of this technology lies in the electrolyte. The electrolyte is a caustic water solution with 25%-30% of potassium hydroxide (KOH). Also can be used as electrolyte, sodium hydroxide (NaOH) or sodium chloride (NaCl). The electrolyte is also used as the catalyst to promote electrolysis to take place. The ions produced in cathode side are transported to the cathode thanks to the liquid electrolyte. The principle of operation is quite simple. The reactions that involve the process are shown in Equation 2and Equation 3. 2H2O+ 2e−→H2+ 2OH−(2) 2OH−→1 2O2+H2O+ 2e−(3) As can be seen in Equation 2, in the cathode electrode water molecules are reduced to hydrogen gas when the current is applied to water. At the same time, as can be seen in Equation 3, at the anode, oxygen arises and a water molecule is generated at the same time. As has been mentioned, the electrolyte has to separate the gases that are produced in each electrode, but this is not possible to be done with a liquid electrolyte. This is why it is necessary to use a separator also known as diaphragm, which has to fulfill the same requirements as a solid electrolyte. In Figure 1, a schematic representation of the explained operation is shown. Figure 1: Schematic representation of the operation of an alkaline electrolysis cell [8]. Thermal management study of an electrolytic stack using CFD page. 17 Some of the technical characteristics of this technology are the following [9]: 1. High efficiency: Is a highly efficient method for producing high-purity hydrogen, with an energy efficiency approximately up to 80%. 2. Low operating voltage: It operates with a relatively low operating voltage of around 1.82.2 volts. This helps to reduce energy consumption and increase efficiency 3. High purity hydrogen: Alkaline hydrogen electrolysis produces high-purity hydrogen gas of up to 99.999%, making it suitable for a wide range of industrial applications. 4. Durability: Alkaline hydrogen electrolysis cells are typically made of durable materials, such as nickel or stainless steel, which can withstand the harsh conditions of electrolysis. However, although this technology has many advantages, it also has some disadvantages, for example [9]: 1. Corrosion: The highly alkaline environment of the electrolyte solution can lead to corrosion of the electrodes and other components of the electrolysis system, which can reduce its lifespan and increase maintenance costs. 2. Limited flexibility: Alkaline hydrogen electrolysis systems are typically optimized for a specific production rate and operating conditions, which can limit their flexibility to adapt to changing demand or energy input. 3. Expensive capital cost: Alkaline hydrogen electrolysis systems can have high capital costs due to the need for high-quality materials and sophisticated control systems to operate the process efficiently. 3.1.2 Proton Exchange Membrane electrolyzer This type of technology is newer than the one mentioned above, but commercial applications of PEM electrolyzers already exist nowadays. The basic principle of operation is the same as for any electrolytic cell, but the main differences lie in the electrolyte used and the configuration of this type of cells. Here, the electrolyte membrane is typically made of perfluorinated sulfonic acid polymer, which conducts protons but is impermeable to gas. So in this technology, the hydrogen protons crosses the electrolyte from the cathode to the anode. The configuration of the electrolytic cell is the following: the electric current reaches the cell through the bipolar plates, an electrically conductive component of the cell. These components are on both the anode and cathode side. In addition, these bipolar plates have channels that page. 18 Report serve to supply water and collect the products produced. After the BPs, the next element to appear is the Porous Transfer Layer (PTL), also on both sides (anode and cathode). Finally, the last layer appearing on both sides is the Catalyst Layer (CL). The reason for having the PTL is to prevent water flowing through the channels from coming into contact with the hydrogen protons that have not yet crossed the electrolyte when the electrolysis reaction is taking place. The configuration of the PEM cell can be seen in Figure 2 Figure 2: Schematic representation of the configuration and operation of a PEM electrolysis cell [10]. For the operation of PEM electrolyzers, water is fed into the anode through the channels located in the BP. When it comes into contact with direct current, this water is splited into oxygen gas and positively charged hydrogen ions (protons) by the oxidation reaction shown in Equation 4. The positively charged protons migrate through the electrolyte membrane to the cathode compartment, where they are combined with electrons from the cathode to form hydrogen gas. The reduction reaction that takes place in the cathode side is shown in Equation 5. 2H2O→O2+ 4H++ 4e−(4) 4H++ 4e−→2H2(5) Some of the technical characteristics of this technology are the following [11]: Thermal management study of an electrolytic stack using CFD page. 19 1. Efficiency: Electrolysis using PEM electrolyzers has a very high conversion of electricity to hydrogen, with an efficiency of around 80-90%. 2. Purity: The hydrogen obtained at the cathode output is very pure. This is due to the fact that the membrane used only allows the hydrogen to pass through, thus preventing any impurities that may be present in the water used from passing to the cathode side. 3. Scalability: Today it is a technology that has already been scaled up commercially. Cells can be connected in series or in parallel to increase production capacity. 4. Rapid response: It can be adapted to the energy demand very quickly because of its high response time. In this way, it can respond quickly to changes in the electric grid. 5. Compact size: PEM electrolyzers are compact, lightweight and modular systems, allowing them to be installed anywhere. This helps the decentralisation of this technology, allowing it to be integrated in places where renewable energy generation technologies are available. 6. Safety: The hydrogen production process using this technology is safe, as it works at low pressures and temperatures. Moreover, it does not require any chemicals that could be harmful. Although these devices have many advantages, the technology is still under development and has the following disadvantages [11]: 1. Cost: Compared to alkaline electrolyzers, this technology is relatively more expensive. This is mainly due to the use of precious materials (Pt, Ir, Ru) which act as electrocatalysts in the electrolysis process. 2. Durability: PEM electrolyzers are devices that require severe maintenance to ensure the durability and stability of these systems. 3. Sensitivity to impurities: This technology is very sensitive to impurities in the water at the cathode inlet. The presence of these impurities can reduce the efficiency of the membrane. This is why this technology requires a pre-treatment process of the water used. 3.1.3 Solid Oxide Electrolyzer This type of technology for producing water electrolysis is the newest technology available today and despite its promising applications is still under development. The most challenging part of the manufacturing of this device is the materials to be used, as materials have to be found page. 20 Report that guarantee the durability of the electrolyzer when working under extreme conditions of corrosion and temperature. Since the present project will be based on this type of electrolyzer, the functions and materials of each layer of the electrolyzer are explained in more detail in Section 4.1. This type of electrolysis is also known as high-temperature electrolysis, as the steam entering the anode has a temperature of between 700ºC and 900°C. By operating at such a high temperature, the idea is that such electrolyzers can work in conjunction with other types of power generation technologies that generate steam at such high temperatures, such as fourth-generation nuclear power plants. Once the electricity has been generated from the high temperature steam generated in the reactor, this remaining steam can be fed to electrolyzer stacks where hydrogen is generated which can then be converted back into electricity at peak demand times. Figure 3 shows a graphical representation of this electrolyzer. Figure 3: Schematic representation of the configuration and operation of a SOEC electrolysis cell [12]. The principle of operation of this technology is as follows: steam enters the inlet of the cathode, which in turn is fed by an electric current to its electrode. At the cathode, the reduction reaction shown in Equation 6takes place. After this reaction, hydrogen is generated which is collected at the outlet. H2O+ 2e−→H2+O2−(6) On the other hand, at the anode side the oxidation of the oxygen takes place. Oxygen ions Thermal management study of an electrolytic stack using CFD page. 21 generated in the cathode, cross the electrolyte and and recombine with the electrons from the anode, causing the oxygen ions to convert back to oxygen. This oxygen is picked up by the air entering through the anode channels, which is expelled through the outlet of the anode channels. The reaction that takes place in the anode is shown in Equation 7. O2−→O2+ 2e−(7) Some of the advantages of this technology are the following: 1. High efficiency: This kind of electrolyzer operates with high temperatures, this means that less electrical energy has to be used to produce the electrolysis. 2. Potential for co-generation: As explained above, this technology can utilise the steam produced in other energy sources such as nuclear or geothermal, allowing the efficiency of these to be increased and electrolyzers to be used as a form of energy storage. 3. Reversibility: SOEC can be operated in both electrolysis mode for hydrogen production and fuel cell mode for power generation, providing a flexible and efficient energy storage solution. However, this technology is still under research and it has some disadvantages. Some of them are the following: 1. High cost: SOEC type electrolyzers use materials that are very expensive, which is why SOEC electrolyzers are more expensive than alkaline electrolysis. 2. High temperature operation: While high operating temperatures allow for high efficiency, they also require more energy to reach and maintain, which can increase overall operating costs. 3. Degradation of materials: The high temperature and corrosive conditions found in this type of electrolyzer mean that the materials used are easily degraded. This is why it is necessary to use materials that can withstand these conditions and guarantee the durability of the device. 4. Emerging technology: This technology is still under development and there are still many barriers to overcome before its commercial application becomes a reality. Its commercial applications are still very limited. page. 22 Report 3.2 OpenFOAM OpenFOAM is a free and open-source CFD software package that allows users to simulate complex fluid flows using numerical methods. In order to obtain the results of this project, OF has been used from a virtual machine. It has been decided to use this and not other software such as Ansys because OF is highly modular, which allows customizing according to the user’s interests. This is why it is widely used in the research field [13]. However, it is not easy to start working with it due to the lack of a graphical interface. OF consists of a collection of solvers and utilities that can be used to model a wide range of fluid dynamics problems, including laminar and turbulent flows, multiphase flows, incompressible and compressible flows, heat transfer and more. Specifically, as will be seen later, the solver dedicated to solving heat transfer will be used to develop this project. 3.2.1 Fundamentals OF is a software that works with a PISO (Pressure Implicit with Splitting of Operators) type algorithm, which plays an essential role in fluid analysis [14]. This algorithm is used in OF to solve the Navier-Stokes equations, as it allows the complexity of the flow simulations to be solved by splitting the temporal and spatial resolution of the equations. The Navier-Stokes equations are fundamental to describe the conservation of mass and momentum in a fluid, and are of great importance in the mathematical formulation of physical problems since they model physically how a flow behaves when forces act on it. These equations are three: conservation of mass, momentum and energy. Equations 8,9and 10 show these conservation laws written in generic form for mass, momentum and energy respectively [15]. ∂ρ ∂t +∇ · (ρu)=0 (8) ρDu Dt =∇ · σ+ρb (9) ρDe Dt +ρDK Dt =−∇ · q+ρr +∇ · (σ·u) + ρb ·u(10) In order to develop these equations and apply them to the project to be developed, it is necessary to apply three concepts. The first of these is the Newtonian fluid, which applies to most fluids. It Thermal management study of an electrolytic stack using CFD page. 23 states that a fluid at rest (or uniform velocity) does not sustain shear stress τ; it can be expressed by Equation 11 [15]. τ= 2µdevD and σ=τ−pI (11) The second concept is the Fourier’s law, which describes heat transfer through a material and states that the rate of thermal conduction is proportional to the temperature gradient in the material. This law is expressed by Equation 12 [15]. q=−k∇T(12) The last concept to be introduced is the material time derivative. This is a is a concept in fluid dynamics that represents the rate of change of a property of a fluid particle as it moves through space, accounting for both spatial and temporal variations. The expression for the material time derivative is shown in Equation 13 [15]. Dψ Dt =∂ψ ∂t +u· ∇ψ(13) Using these three concepts and and developing the Navier-Stokes equations for a compressible, Newtonian fluid (which is the case in this project), the following equations are obtained: conservation of mass, momentum and energy (Equations 14,15, and 16). ∂ρ ∂t +∇ · (ρυ) = 0 (14) ∂ρυ ∂t +∇ · (ρυυ) = −∇ρ+∇ · ¯τ+ρg (15) ∂ρe ∂t +∇ · (ρυe) = ∇ · (µ∇T)−∇·(pυ) + ∇ · (¯τ·υ) + ρg ·υ (16) Finally, Equation 17 shows the equation of state of ideal gases, which is also needed for fluids problem solving and models the behaviour of many gases under typical working. page. 24 Report pV =nRT (17) It is also important as fundamentals of the CFD discipline to explain the general transport equation, which describes how a physical property is transported in a spatial domain over time. This is shown written in the generic form in Equation 18. ∂ϕ ∂t +∇ · (vϕ) = ∇ · (k∇ϕ) + S(18) The first term of this equation makes reference to the partial derivative of the property ϕwith respect to time. The second term is known as the advective term, and refers to the transport of properties of a fluid due to its movement. The third term of the Equation 18 is the diffusive term, and makes reference to the process by which particles are mixed and uniformly distributed in a fluid due to random molecular motions. Finally, the last term is the source term, and represents sources or sinks of property ϕin the domain. 3.2.2 Case structure As it has been explained, OF is a very powerful tool to solve multiple physical problems using CFD simulations, but one of its main disadvantages when starting to use it is that it lacks a graphical interface, which makes it difficult to understand how proceed to start a simulation. The aim of this section is to explain the basic concepts of OF in order to understand how the simulation of the project has been carried out. This will allow: •Understand the underlying folder structure of each OF simulation within the context of this project. •Set up the key files and accurately verify the completion of each aspect, thus ensuring the successful launch of the simulation in the context of this project. To start the process, it is suggested to start from a previously solved case, which will be copied to a new folder. An effective approach involves identifying the type of phenomenon to be simulated. This solved case can be found as a tutorial set by OF in the installation directory. Then, it can be proceed to duplicate the relevant file. It is advisable to assign the name "base case" to this instance in order to distinguish it from future studies that may be derived from it. The general structure of an OF case is the following: •0 folder: This folder stores the initial conditions for each simulation. Several files are Thermal management study of an electrolytic stack using CFD page. 31 Mainly, as can bee seen in Figure 7, a SOEC is typically constituted of four layers: The hydrogen or also named as steam electrode (which is the cathode), the electrolyte, the oxygen or air electrode electrode (which is the anode), the intermediate layer (also called barrier layer), which is located between the electrolyte and the oxygen electrode, and the interconnector. This last layer is more know as Bipolar Plate (BP). When there is a stack of cells (more than 2 cells), the interconnector allows to joint the two cells. 4.1.1 The hydrogen electrode (Cathode) In the hydrogen electrode is where the steam is electrochemically reduced into hydrogen. This layer can be divided in two sub-layers, the structural and the functional part. The structural part has a more constructive function, and is responsible for joining the bipolar plate with the functional part of the hydrogen electrode. Hydrogen reduction can also occur in the structural part of the hydrogen electrode, but it differs from the functional part in that its porosity is greater. This means that when the reduction occurs and the molecules are divided into hydrogen and oxygen, the process can be reversed when they come into contact with other water molecules. On the other hand, the functional part of the hydrogen electrode is less porous, to prevent the reaction from being reversed. It is in this part that the water reduction reaction is assumed to occur. The reaction takes place at the interface of the hydrogen electrode and the electrolyte, also known as the Triple Boundary Layer. It is named as such because it is where the electrical phase (electron conductive), the ionic phase (O2−conductive) and the gaseous phase (steam supply and hydrogen release) meet [1]. As properties to highlight, this layer has to allow the conduction of electrons and ions. Also, it must have catalytic properties so that the reaction can take place, and allow the transport of water, hydrogen and oxygen. In addition, it has to withstand the high temperatures of the steam being used. Currently a metal-ceramic composite of the material used in the electrolyte (yttriastabilized zirconia) and Nickel is used. This material is a non-precious metal catalyst with high electronic conductivity. In Equations 19,20,21 and 22, the density (ρ), specific heat capacity (Cp), thermal conductivity (k) and electrical resistivity (designated as αin this document) in function of the temperature are shown for YSZ-Ni. Also in Figures 8,9,10 and 11 a graphical representation of these properties are shown. The properties of this material have been obtained from a research paper where the thermo-physical properties [19] of this material were analyzed. To plot the figures of this document with enough resolution and obtaining the equations to calculate each parameter in function of temperature, the data from the original graph (the one that appears in the research paper) have been extracted making use of the software tool WebPlot- page. 32 Report Digitizer. This data was then entered into Excel to draw the graph and obtain the polynomial that best fits it, thus obtaining the equation to represent each property as a function of temperature. In these graphs the dotted line is the mathematical fit to the experimental values found in the article. This procedure is the one that has been carried out for all other materials in this section (for all other layers of the electrolyzer). ρ(T) = −8·10−15 ·T5+ 4 ·10−11 ·T4−6·10−8·T3+ 4 ·10−5·T2−0.008 ·T+ 4.8665 (19) Cp(T) = −8·10−15 ·T5−3·10−11 ·T4+ 4 ·10−8·T3−3·10−5·T2+ 0.0074 ·T−0,2031 (20) k(T) = 1 ·10−13 ·T5−3·10−10 ·T4+ 3 ·10−7·T3−0.0001 ·T2+ 0.0143 ·T+ 5.8851 (21) α(T) = 38.559 exp−0.004T(22) Figure 8: Density in function of temperature for YSZ-Ni. Thermal management study of an electrolytic stack using CFD page. 33 Figure 9: Specific heat in function of temperature for YSZ-Ni. Figure 10: Thermal conductivity in function of temperature for YSZ-Ni. Figure 11: Electrical resistivity in function of temperature for YSZ-Ni. page. 34 Report 4.1.2 Electrolyte The electrolyte is the layer of the SOEC that determines the composition of the surrounding layers (hydrogen electrode, oxygen electrode and intermediate layer). This layer has to avoid the ionic conduction between the cathode and the anode, concretely oxygen ion conduction in the case of a SOEC electrolyzer. At the same time, the layer has to avoid the conduction of electrons (electrically insulating) and the transport of gasses between the anode and the cathode (needs to be enough dense), in this case, the transport of the steam. Moreover, the electrolyte has to sustain the extreme environment that evolves the process of electrolysis in this kind of electrolyzers. This means that the material has to have good mechanical strength, chemical and thermal resistance properties in order to ensure the durability of the device. Currently, yttria-stabilized zirconia (YSZ) is the material used as electrolyte, as it is the one that has shown good and stable performance over time in the typical SOEC temperature range 700–850 ºC thanks to a high ionic conductivity associated with good thermal and chemical stabilities [12]. In Equations 23,24,25,26, the density (ρ), specific heat capacity (Cp), thermal conductivity (k) and electrical resistivity (α) in function of the temperature are shown for YSZ. Also in figures 12,13,14 and 15, a graphical representation of these properties are shown. The properties of this material have been obtained from research papers where the thermo-physical properties [19] [20] of this material were analyzed. To plot the figures of this document with enough resolution and obtaining the equations to calculate each parameter in function of temperature, the data from the original graph have been extracted making use of the software tool WebPlotDigitizer. ρ(T)=2·10−14 ·T5−6·10−11 ·T4+ 6 ·10−8·T3−2·10−5·T2+ 0.0016 ·T+ 5.5644 (23) Cp(T) = −2·10−7·T2+ 0.0004 ·T+ 0.4706 (24) k(T) = −1·10−7·T2+ 0.0002 ·T+ 1.9748 (25) α(T) = +0.0202 ·T2−41.376 ·T+ 21167 (26) Thermal management study of an electrolytic stack using CFD page. 35 Figure 12: Density in function of temperature for YSZ. Figure 13: Specific heat in function of temperature for YSZ. Figure 14: Thermal conductivity in function of temperature for YSZ. page. 36 Report Figure 15: Electrical resistivity in function of temperature for YSZ. 4.1.3 The intermediate layer The intermediate layer connects the electrolyte with the oxygen electrode, and is used to avoid thermal mismatching between these two layers. This thermal mismatching can occur due to the materials of which the oxygen electrode is normally made, since high temperatures can produce an expansion of the thermal coefficient, in the case of electrodes made of cobalt, thus producing an increase of the ohmic resistance. Due to this fact, this intermediate layer can be understood as a transition layer to allow good thermal compatibility and avoid chemical interactions (element migration) while ensuring a good ionic conductivity [12]. In addition to the characteristics mentioned above, this layer must have catalytic properties and allow the transport of a gas phase. Gadolinium-doped cerium oxide (CGO or GDC) and yttrium-doped cerium (YDC) are common reference materials. In this work, to obtain the properties, it will be assumed that the layer is composed of GDC. In Equations 27,28,29 and 30 the density (ρ), specific heat capacity (Cp), thermal conductivity (k) and electrical resistivity (designated as αin this document) in function of the temperature are shown for GDC respectively. Also in Figures 16,17,18 and 19 a graphical representation of these properties are shown. The properties of this material have been obtained from a research paper where the thermo-physical properties of cerium alloys were studied [21]. To plot the figures of this document with enough resolution and obtaining the equations to calculate each parameter in function of temperature, the data from the original graph have been extracted making use of the software tool WebPlotDigitizer. ρ(T) = −3·10−6·T2+ 0.0072 ·T+ 0.8905 (27) Thermal management study of an electrolytic stack using CFD page. 37 Cp(T) = 6 ·10−15 ·T5−1·10−11 ·T4+ 7 ·10−9·T3−2·10−6·T2+ 0.0003 ·T−0.3495 (28) k(T)=7·10−9·T3−7·10−6·T2+ 0.0034 ·T+ 5.3779 (29) α(T) = −5·10−6·T4+ 0.0086 ·T3−5.2539 ·T2+ 1373.1·T−125734 (30) Figure 16: Density in function of temperature for GDC. Figure 17: Specific heat in function of temperature for GDC. page. 38 Report Figure 18: Thermal conductivity in function of temperature for GDC. Figure 19: Electrical resistivity in function of temperature for GDC. 4.1.4 The oxygen electrode (Cathode) In the oxygen electrode is where the oxygen ions (O2−) are oxidized to oxygen. Normally, this reaction takes place in the interface of the electrolyte with this electrode. This oxygen produced is collected by the air entering the inlet of the cathode channels. The main problem with the connection between these two layers is the thermal and chemical compatibilities with the materials used. Here remains the importance of the intermediate layer, explained in the section before. In addition to the fact that the material used as the oxygen electrode must allow for such compatibility, other required properties are that it is a good ionic and electrical conductor. It must also be a good catalyst to encourage the oxygen oxidation reaction to take place. Also has to Thermal management study of an electrolytic stack using CFD page. 39 allow the transport of gasses. Currently the materials that are used are cobalt or strontium. For this project it will be assumed that the layer is made of Cobalt. In Equations 31,32,33 and 34 the density, specific heat, thermal conductivity and electrical resistivity in function of the temperature are shown for Cobalt respectively. Also in Figures 20,21 and 22 a graphical representation of the specific heat capacity, thermal conductivity and electrical resistivity in function of the temperature is shown. The properties of cobalt have been obtained from a research paper [22].In this case, no graphical representation of the density is shown as the formula was obtained directly from a scientific article [23]. ρ(T) = 1 1.1997 ·10−5·T+ 0.10723 (31) Cp(T)=0.0002 ·T+ 0.421 (32) k(T)=3·10−13 ·T5−7·10−10 ·T4+ 6 ·10−7·T3−0.0002 ·T2+ 0.0261 ·T+ 13.652 (33) α(T) = −8·10−19 ·T4+ 2 ·10−15 ·T3−2·10−12 ·T2+ 8 ·10−10 ·T−7·10−8(34) Figure 20: Specific heat in function of temperature for Cobalt [22]. page. 40 Report Figure 21: Thermal conductivity in function of temperature for Cobalt [22]. Figure 22: Electrical resistivity in function of temperature for Cobalt [22]. 4.1.5 The interconnector The interconnector, also known as Bipolar Plates, plays an essential role when more than two electrolytic cells are joined together (a stack). It is what allows one cell to be joined to the other, thus separating the oxygen electrode of the first cell from the hydrogen electrode of the second cell. In addition, the electric current from one cell to its neighbouring cell is transferred through this layer. In addition, it is through the channels in the BPs that the inlet steam is fed to the cathode and the inlet air to the anode. The material from which the BPs have to be made has to be able to withstand the high inlet steam temperatures (700-800°C) and be resistant to the corrosive conditions that can occur during electrolysis. On the other hand, it must be a good conductor of electrical current and have a Thermal management study of an electrolytic stack using CFD page. 47 ¯ k=     ¯ kx0 0 0¯ ky0 0 0 ¯ kz      (40) To calculate the thermal conductivity in each direction, so-called equivalent thermal resistance circuits are used. The thermal resistances for the case study are shown in Figure 27. Looking at this figure, it can be seen that the equivalent thermal resistance for the x and y direction is the sum of their parallel resistances, while for the z direction, it is the sum of their series resistances. The generic expression for both sums (sum in parallel and sum in series) is shown in Equations 41 and 42 respectively. In this equations, ∆zirefers to the thickness of the layer, kiis the thermal conductivity of the material and Ai is the cross sectional area perpendicular to the direction of the conduction in the layer. 1 Rx,y =¯ kx,yAx,y ∆z= n X i=1 1 Ri = n X i=1 kiAi ∆zi (41) Rz=∆z ¯ kzAz = n X i=1 Ri= n X i=1 ∆zi kiAi (42) By substituting terms and rearranging the Equation 41 for the repeating element illustrated in Figure 27, one can calculate the equivalent thermal conductivity in both the x-axis and y-axis. Repeating the same but with Equation 42, the equivalent thermal conductivity in z-axis can also be calculated. The resulting expressions are shown in Equations 43,44 and 45 for directions x, y and z respectively. ¯ kx=2ksteel∆zIC +kNi,f ∆zNi,f +kcell∆zcell +kCu,f ∆zCu,f ∆zRE (43) ¯ ky= 2ksteel ∆IC −∆zch wch wRE +kNi,f ∆zNi,f +kcell∆zcell +kCu,f ∆zCu,f ∆zRE (44) ¯ kz=∆zRE 2∆zch ksteel wRE (wRE −wch)+ 2∆zIC ksteel +∆zNi,f kNi,f +∆zcell kcell +∆zCu,f kCu,f (45) Going back to equation 39, the first term of this equation refers to the advective part and the page. 48 Report second term to the diffusive part (where the average thermal conductivity calculated previously will be used). This is only applicable in CFD simulations when a fluid is simulated. However, as explained previously, no fluid was simulated in this reference article. Therefore, in order to take into account the effect of the fluid, the researcher imposed a boundary condition known as convection. Convection can be understood in CFD as a boundary condition that involves two phenomena: on the one hand, the transfer of energy associated with random molecular motion, known as diffusion, and on the other hand, the energy transferred through the global or macroscopic motion of the fluid, known as advection. Convective heat transfer is predetermined by Newton’s cooling law, which is shown in Equation 46. dQ dt =hA(Ts−T∞)(46) In this equation, the term T∞refers to the temperature of the fluid, Tsis the temperature of the surface of the solid body, Ais the heat transfer area and hrefers to the convection heat transfer coefficient, which depends on the fluid properties, the fluid velocity and the surface characteristics. Since no fluid is modelled in the theoretical model, convection will be defined as a volumetric source term Sc T, and describes the amount of heat exchanged thorough the surfaces, named as Ach, in the volume of a repeating element (Ω). In Equation 47 is shown the expression used in the reference article to calculate the term Sc T, where his known as the convective heat transfer coefficient and Ach is the area of the solid domain being wet by a specific fluid (the channel’s surfaces areas). This is applicable for a stationary media. Sc T=hAch |Ω|(Tsolids −Tgas)(47) Variable h(heat transfer coefficient) of Equation 47 can be obtained using the Nusselt number. The expression to obtain the Nusselt number is shown in Equation 48, where the kgas makes reference to the thermal conductivity of the fluid. This calculation of the variable hallows the Equation 47 to be rewritten and converted into the Equation 49. In the treated SOEC stack, where laminar flow is fully developed within a rectangular channel, the Nusselt number attains a constant value. Nu =hDh kgas (48) Thermal management study of an electrolytic stack using CFD page. 49 Sc T=Nukgas Dh Ach |Ω|(Tsolids −Tgas)(49) On the other hand, the heat produced by the electrolysis reaction and various losses from the electrochemical reactions and ohmnic losses due to charge transfer is accounted with the heat source term Sr T. This source can be calculated with the Equation 50. This heat is uniformly generated at each point within the model, distributed evenly across the height of each repeating element ∆zRE. In Equation 50,∆Srefers to the entropy change for the reaction, Fis the Faraday’s constant, izis the current density and ASR is the area specific resistance and it will be described in the charge transport equations. Sr T=1 ∆zRE izT 2F+i2 zASR=izT ∆zRE2F+i2 zASR ∆zRE =Ser T+SJoule T(50) The amount of heat generated or consumed is determined by the balance of the two terms in Equation 50. In steam electrolysis, where the reaction heat is negative, there comes a point when it balances with the contributions from the second term, specifically during thermo-neutral operation. It is important to note that the heat source term, Sr T, is expressed per volume. Therefore, it is necessary to divide the current density by the height of the repeating element to convert it into a volumetric measure. It is also important to point out that the first term of this equation (Ser T) refers to the heat generated by the electrolysis reaction, which is only present in the MEA region, while the second term (SJoule T) refers to the heat generated by the Joule effect, which takes place in all regions. Summarizing, the three heat source terms accounting for both convective heat transfer, Sc T, and reactions, Sr tcan be applied for each region (cathode bipolar plate, MEA and anode bipolar plates). The total source term for each region is shown in Equations 52,53,54. ST=Sr T+Sc T(51) SBP,cathode T=hAch |Ω|(Tsolids −Tfuel) + iz ∆zRE (izASR)(52) SBP,anode T=hAch |Ω|(Tsolids −Tair) + iz ∆zRE (izASR)(53) page. 50 Report SMEA T=−hAch |Ω|(Tsolids −Tfuel)−hAch |Ω|(Tsolids −Tair)+iz ∆zRE T∆S 2F+izASR(54) Electrical charge transport To begin with the equations representing the electrical charge transport in a SOEC, it is necessary to start with the equation representing the voltage difference, ∆Vover a repeating element of height ∆zRE. This is expressed in Equation 55, where the term izrefers to current in the z direction, Eis the cell open circuit voltage and ASR represents the area specific resistance. This last term includes the different kinds of resistances that can appear in the repeating element (diffusion, charge transfer, ohmnic losses, etc.). ∆V=E−izASR (55) Following this equation, the stack voltage can be obtained by adding the voltage differences for each repeating element in the stack height. Knowing that the average change in potential over the height is ∆V/∆zthe gradient in the z direction becomes ∂V/∂z. Applying this to Equation 55 and rearranging it to obtain an expression for the variable iz, gives the Equation 56. iz=E ASR −∆zRE ASR ∂V ∂z (56) On the other hand, the expression to obtain the cell open circuit voltage, E, is defined in Equation 57. In this equation, the terms πeq air and πeq fuel refers to the electromotive potentials at equilibrium for each electrode, π∗ iis the chemical potential of species i giving a reference state, the term aiis the activity of the component i in the corresponding gas mixture, ∆Grx,H2Orefers to Gibbs energy change for the water splitting reaction, p−is the standard pressure (1 atm) and finally, pirefers to the partial pressure of the component i in the corresponding gas mixture. E=πeq air −πeq fuel = 1 2µ∗ O2+µ∗ H2−µ∗ H2O 2F+RT 2Fln    aH2a 1 2 O2 aH2O  (57) There are different expressions for the term ASR. For this explanation, it has been decided to use the model of Leonide et al [24], which expresses the ASR as a function of conecntration overpotential, ηconc, and the activation overpotential, ηact, for both electrodes, and also in function Thermal management study of an electrolytic stack using CFD page. 51 of the Ohmnic overpotential, ηOhm. With this expression, the cell voltage ∆V, can be expressed as it is shown in Equation 58. In this equation, f and o denote the fuel/oxygen electrode respectively. ∆V=E−(ηOhm +ηact,f +ηact,o +ηconc,f +ηconc,o)(58) Expressions to obtain the different overpotentials mentioned before can be obtained in reference [24]. This expressions does not take into account the resistance of the interconnectors, coating, in-plane conduction of current through the contact layers for current collection, etc. For this reason, 0.2 Ωcm−2has been added to this ASR expression in the modelling [5]. Following the explanation, one of the laws that must be fulfilled when modelling charge transport is that the current must be conserved over the entire computational domain. This is expressed in Equation 59, where irepresents the current generated by the electrical field and the external current density, also known as the current density flux. ∇ · i= 0 (59) In the direction perpendicular to the plane (z-direction), Equation 59 can be transformed into Equation 60. In this equation, the term E/ASR can be identified as the external current density and ∆zRE refers to the conductivity in the z-direction. Within the cell plane, specifically in the xand y-directions, no "external current" is generated, and the conductivity, denoted as γ, is the average of the unit cell. This average conductivity is primarily influenced by the interconnect, which is considerably higher than the out-of-plane conductivity. Consequently, Equation 59 transforms into Equation 61 for the x-direction and Equation 62 for the y-direction. goo ∂ ∂z E ASR −∆zRE ASR ∂V ∂z = 0 (60) ∂ ∂x −γx ∂V ∂x = 0 (61) ∂ ∂y −γy ∂V ∂y = 0 (62) In relation to the above equations, the conductivity tensor and the external current density vec- page. 52 Report tor can be expressed as shown in Equations 63 and 64 respectively. Finally, the charge transfer equation can be rewritten as Equation 65 γ=     γx0 0 0γy0 0 0 ∆zRE ASR      (63) ie=     0 0 E ASR      (64) ∆·(γ∆V)=∆·ie(65) 4.3.2 Simplified theoretical model equations The equations explained in the previous section are too complex to be replicated in the OF model to be developed, as it is beyond the scope of this project to develop such a detailed model. For this reason, in order to simulate the behaviour of a SOEC when it is in operation, some of the equations that have been explained have been taken and some hypotheses and simplifications have been made in order to apply them. It is important to note that OF has solvers that allow to solve directly some of the equations that have been set out above. For this project, the solver to be used is the one known as chtMultiRegionSimpleFoam (CMRSF). CMRSF is a solver that allows to solve conjugate heat transfer problems between multiple regions, which allows to simulate the coupling between a fluid flow and a solid, and therefore to solve heat transfer in different regions and domains where one of the regions is a solid and the other a fluid. Which equations the solver will solve depends on the case configuration. That is, depending on which regions are declared as solid and which others as fluid. Other important simulation parameters to set at the beginning are the initial and boundary conditions in 0/, and the definition of the material properties in constant. Returning to the equations to be solved, the equations to be implemented to model the heat transfer in a SOEC are the Equations 53,52 and 54. In these equations, the most complex part to implement in OF is the last term, which represents the heat generated by the electrolysis reaction, Ser T, and the heat generated by the Joule effect SOhm T, since OF does not have these Thermal management study of an electrolytic stack using CFD page. 53 equations in the solver to be used. These two terms are represented separately in the Equation 66. Ser T=izT∆S ∆zRE2F, SJoule T= i2 zASR ∆zRE !(66) As it has been done in the explanation of the theoretical equations (Section 4.3.1, for this part it will be explained how the fluid transport, heat transfer and electrical charge transport are modelled in the project. The heat transfer part is where the implementation of the term Ser Tin OF has been addressed. On the other hand, in the electrical charge transport part, it is explained how the Joule effect losses, to which the term SOhm Trefers, have been introduced in OF. Fluid transport Some simplifications will be taken to model fluid transport. On the one hand, as mentioned above, it will be assumed that no part of the electrolyzer is porous, so the fluid will only interact with solid media. Furthermore, the influence of the fluid on species transport will be simplified by assuming perfect (ideal) transport, i.e. all the produced oxygen reaches the oxygen channel and no reverse flow or crossover hydrogen/water is assumed. It will also assume an ideal electrolyzer operation, where all inlet steam is converted to hydrogen and no steam remains in the outlet. Finally, it is necessary to point out that when using OF and indicating which regions are a fluid, the program itself will solve the corresponding mass and momentum conservation equations for each of the fluids at the cathode and anode. In other words, the equations explained in Section 3.2.1 will be solved for the fuel/oxygen in the computational domain defined. Heat transfer As explained in Section 3.2.1, one of the advantages of OF is that it allows to solve the NavierStokes equations. This is one of the advantages of CFD, since it allows solving the advective and diffusive term of Equation 39 without the need to model a boundary condition such as convection (as was done in the theoretical model), since the fluid is simulated. All that has to be done is specify which regions are solid and which regions are fluid. Apart from the common heat transfer mechanisms explained before, it is also necessary to include in the model the equation representing the heat that is generated when electrolysis takes place (source term Ser Tin Equation 66). This equation is not included in the OF solvers and therefore it is necessary to find a way to introduce it into the model. After investigating dif- page. 54 Report ferent ways of doing this, it has been decided to introduce it into the model developed in this project as a source term within the region representing the membrane (with all its layers included within the region). Using the fvOptions option within the 0/T file within the corresponding region, a constant heat source can be simulated. The value of such a heat source shall be calculated using the equation of the source term Ser Tfor the boundary conditions set in the model. That is, this heat source will not change throughout the simulation and will be constant during all the time. The variables T and ∆Sshall be kept constant throughout the simulation, using the value of ∆Sassociated with the boundary condition temperature. This means that if the boundary condition temperature is 600°C, the ∆Svalue to be used is the one corresponding to that temperature. It is important to see that a value of izremains to be added in order to calculate the constant value of the heat source Ser T. To obtain this value of iz, which once calculated will remain constant, a voltage, V, has been imposed. With the resistivity of each layer obtained in Section 4.1, the resistance of each of these layers can be obtained applying Equation 67. Combining the resistances of each layer (in series) it can be obtain an equivalent resistance. Finally, with the equivalent resistance obtained and the imposed voltage, Ohmn’s law (Equation 69) is applied to obtain the value of the current, which will be the iz value to be applied in Equation 66 to obtain the value of the constant heat generated by the electrolysis reaction, Ser T. To illustrate this last explanation, an example will be used. For this example, two layers will be used, layer 1 and layer 2. In the following equations, the steps explained previously can be followed, being Ri,αi,liand Aithe electrical resistance, resistivity, length and transversal area of each layer respectively. Ri =αi li Ai (67) Requivalent = n X i=1 =R1 + R2(68) V=izRequivalent (69) Electrical charge transport Another source that needs to be modelled in this project is the heat dissipated due to the Joule effect. This source is the one that has been defined in Equation 66 as SJoule T. The Joule effect Thermal management study of an electrolytic stack using CFD page. 55 occurs when electric current flows through a material, due to the electrical resistance of the material, and takes place in all the regions defined in the computational domain. As a consequence, energy is lost in the form of heat. This dissipated power can be calculated using equation 70, where Iis the electric current flowing through the conductor and Rrepresents the electrical resistance of the conductor. P=I2·R(70) This equation can also be written as a function of the electric field, E, and the electric current density, J. Applying Equation 71 to obtain an expression for the electric field and Equation 72 for the electric current density and substituting these equations into Joule’s law, it can be obtained Equation 73, which expresses the power dissipated by the Joule effect (expressed previously in Equation 70) as a function of the electric potential (V e) and the electrical conductivity, σ.  E=−∇ · V e (71) J=−σ(T) E=−σ(T)∇V e (72) P=I2·R=J·E= (σ∇V e)· ∇V e (73) There is an option in OF that allows to implement Equation 73 in the simulation, thus representing the contribution of the Joule Heating effect to a thermal solver. The option solves an equation for the electrical potential, V, of the form of the Equation 65, whereas at the same time adds the heat source of the Equation 73. To apply this option, a voltage has to be defined as a boundary condition in the 0 folder. The structure of the JouleHeatingSource option within OF is as shown below. heating { type jouleHeatingSource ; active true ; jouleHeatingSourceCoeffs { page. 56 Report anisotropicElectricalConductivity no; // O pt i o na l l y s p e c i f y sigma as a f u n c t i o n o f t em p er a tu r e sigma 84033613.45; // // sigma t a b l e //( // (0 127599.8469) // (1000 127599.8469) //) ; } } The electrical conductivity of the previous equations can be defined in the option JouleHeatingSource in different ways: •If not present the sigma field will be read from file (standard field sigma). •If the sigma entry is present the electrical conductivity is specified as a (potentially uniform) function of temperature. For example, using a temperature sigma table like the following: ((0 127599.8469)(1000 127599.8469)). In this table, a sigma value is specified for a given temperature. •If the anisotropicElectricalConductivity flag is set to true, sigma should be specified as a vector quantity. Thermal management study of an electrolytic stack using CFD page. 63 { // Surface −wise min and max r ef i n e m e nt l e v e l le v el (2 2) ; } } resolveFeatureAngle 30; planarAngle 30; refinementRegions{} // Mesh s e l e c t i o n // ~~~~~~~~~~~~~~ allowFreeStandingZoneFaces true ; locationsInMesh ( ((5 10 −2) BP1) // c e l l Z o n e 0 ((5 10 5) MEA) // c e l l Z o n e 1 ((5 10 12) BP2) // c e l l Z o n e 2 ) ; faceZoneControls { } } Once the above file has been prepared, it is only necessary to run the command snappyHexMesh -overwrite in the terminal, and OF will start meshing the geometry. Figure 33 shows the result obtained after meshing the BM. Figure 33: Mesh generated with snappyHexMesh for the BM. page. 64 Report The last step is to create the regions. To do this, OF has a built-in tool that, after having previously declared the regions, the software itself will separate them and create the patches automatically for each one of the regions. This tool is called splitMeshRegions and to run it, it is needed to type the command splitMeshRegions -cellZoneOnly -overwrite in the terminal. It is also necessary to have previously defined in the regionProperties file (inside the constant folder) what type of region each one is, whether solid or fluid. Figure 34 shows the result obtained for the BM after applying the splitMeshRegions. It is also important to note that by executing this command, OF creates in each of the main folders (0, constant and system) a folder for each region. That is, in this case, a folder, in each of the main folders, will be created with the name BP1, another with B2 and another with MEA. Figure 34: Regions created with splitMeshRegions for the BM. 5.3 Preparation of the BM simulation The first step before starting to define other files is to have, as aforementioned, prepared the file defining whether a region is a solid or a fluid. This is very important because OF will look for some data or others when running the simulation based on this file. The regionProperties file used for the BM is shown below. FoamFile { version 2 . 0 ; format a s c i i ; cl as s dictionary ; obje c t regionProperties ; Thermal management study of an electrolytic stack using CFD page. 65 } // ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗// regions ( flui d () soli d (BP1 BP2 MEA) ) ; Once the mesh is ready, the simulation files for each region will be prepared. The first file to modify is the thermophysicalPoroperties file, which is located in the constant folder. Table 3shows the properties defined for each region and their corresponding values. Table 3: BM thermophysical properties for each region. Thermophysical property BP MEA Units Description molWeight 55.85 349.03 u.m.a Molecular weight kappa 23.94 2.07 W/mK Thermal conductivity Cp 547.2 652.6 J/kg K Specific heat capacity ρ7546.3 4800 kg/m3Density The properties of region BP1 and BP2 are the same, since they are made of the same material. This is the reason why in the previous table appears as region BP. These properties have been obtained from the data presented in Section 4.1 for a temperature of 700ºC. These properties could be set as a function of temperature, and OF could recalculate them for each temperature value obtained in the simulation. However, to simplify the simulation, they will be kept constant with their initial value for the whole simulation. An example of the thermophysicalPoroperties file for the BP1 region is shown below: FoamFile { version 2 . 0 ; format a s c i i ; cl as s dictionary ; object thermophysicalProperties ; } // ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗// thermoType { type heSolidThermo ; mixture pureMixture ; page. 66 Report transport constIso ; thermo hConst ; equationOfState rhoConst ; specie specie ; energy sensibleEnthalpy ; } mixture { spe ci e {molWeight 5 5 . 8 5 ; } tr ansport {kappa 2 3 . 9 4 ; } thermodynamics { Hf 0 ; Cp 5 4 7 . 2 ; } equationOfState { rho 75 4 6 . 3; } } On the other hand, in the boundary file of each region (inside the polyMesh folder, which is where the mesh is stored), it will be necessary to modify the boundary type for the boundaries between two regions (for expample, BP1_to_MEA). This is done to make the coupling between regions, since when the mesh is created, the coupling between regions is not yet defined. That is, what is being done is to indicate to OF which regions are coupled, so that it is able to read variables between regions. An example of the boundary file for BP1 region with the corresponding modification already done is shown below. As it could be seen, the boundary BP1_to_MEA is defined as a mappedWall type. Here is where the previous explanation is applied. FoamFile { version 2 . 0 ; format a s c i i ; arch "LSB ; l a b el =32; s ca l a r =64" ; cl as s polyBoundaryMesh ; l oc at i on " constant /BP1/polyMesh " ; obje c t boundary ; } // ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗// ( i n l e t { Thermal management study of an electrolytic stack using CFD page. 67 type patch ; nFaces 384; s t a rtFac e 50796; } walls { type wall ; inGroups 1( wall ) ; nFaces 360; startFace 51180; } bottomwall { type wall ; inGroups 1( wall ) ; nFaces 480; startFace 51540; } topwall { type wall ; inGroups 1( wall ) ; nFaces 480; startFace 52020; } BP1_to_MEA { type mappedWall ; inGroups 1( wall ) ; nFaces 6144; start F a ce 52500; sampleMode nearestPatchFace ; sampleRegion MEA; samplePatch MEA_to_BP1 ; } ) It will also be necessary to define the boundary conditions in each region. These boundary conditions are the ones defined above (Section 5.1). The temperature and voltage boundary conditions file for the BP1 region is shown below as an example. As can be seen, the voltage at the input and output of the BP1 region is the same. This is because the electrical resistance of the bipolar plate is considered negligible, as its value is very low. Because of this, there is no voltage drop between these two points. FoamFile { version 2 . 0 ; format a s c i i ; arch "LSB ; l a b el =32; s ca l a r =64" ; cl as s volScalarField ; location " 0/BP1" ; obj ec t T ; } // ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗// dimensions [ 0 0 0 1 0 0 0 ] ; i n t e rn a l F i e l d uniform 2 9 6 . 9 ; page. 68 Report boundaryField { #includeEtc " caseDicts / setConstraintTypes " i n l e t { type fixedValue ; value uniform 700;} walls { type zeroGradient ; } bottomwall { type zeroGradient ; } topwall { type zeroGradient ; } BP1_to_MEA { type compressible :: turbulentTemperatureCoupledBaffleMixed ; value $internalField ; Tnbr T ; kappaMethod solidThermo ; } } FoamFile { version 2 . 0 ; format a s c i i ; cl as s volScalarField ; object jouleHeatingSource :V; } // ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗// dimensions [ 1 2 −3 0 0 −1 0 ] ; in te rn al F i e l d uniform 0 ; boundaryField { i n l e t { type fixedValue ; value uniform 1 . 3 ; } walls { type zeroGradient ; } bottomwall { type zeroGradient ; } topwall { type zeroGradient ; } BP1_to_MEA { type fixedValue ; value uniform 1 . 3 ; } } Finally, the source terms for each region remain to be defined. The first of these will be the source term Ser T, which will be defined as a constant heat source in the MEA region. This term will be calculated as explained previously using the Equation 66 but first a value for izneeds Thermal management study of an electrolytic stack using CFD page. 69 to be estimated. This will be done by applying Ohm’s law, knowing that the total cell voltage is 1.3V and obtaining an equivalent resistance, Req, for the whole cell. Knowing the electrical resistivities for the BP (stainless steel) and for the MEA (assuming all of it is electrolyte material) the electrical resistance for each material is calculated by applying Equation 67. RBP = 1.19 ·10−8Ωm0.1m 6·10−3m2= 1.98 ·10−7≈0 RMEA = 0.7Ωm0.1m 6·10−3m2= 11.6Ω Req =  RBP 1+  RBP 2+RMEA =RMEA = 11.6Ω As mentioned above, the resistance of the bipolar plate to the current is practically zero, so it can be considered negligible. Finally, applying Ohm’s law, the value of Iis obtained. This value obtained must be expressed in A/m2, so the value obtained is divided by the area through which the current flows (A= 9 ·10−3m2). V=I·Rqu →1.3V=I·11.6Ω →I=1.3V 11.2Ω = 0.11A iz=I A=0.11 9·10−3m2= 12.22A/m2 With this data and knowing that the value of T·∆Sfor 700°C is 55kJ/mol, the Faraday constant, F, is 96485C/mol and the value of ∆RE is 0.2m, the value of the source term due to the electrolysis reaction can be obtained as shown below. Ser T=izT∆S ∆RE2F=12.22A/m2·55000J/mol 96485C/mol ·0.2·2= 17.42W/m3 Moreover, the source term due to the Joule effect applies in all regions. However, it will be relevant in the MEA region, since it is the only one that opposes resistance to the electric current. Both source terms are indicated in the fvOptions file (inside the system folder) for each region. The fvOptions file for the MEA region is shown below. In this file, for the scalarSemiImplicitSource the selected volumeMode is specific, which allows to express a volumetric source. page. 70 Report FFoamFile { version 2 . 0 ; format a s c i i ; cl as s dictionary ; location " system " ; obje c t fvOptions ; } // ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗// MEA_AbsoluteEnergySource { type scalarSemiImplicitSource ; active true ; selectionMode cellZone ; cellZone MEA; volumeMode s p e c i f i c ; sources {h ( 17.42 0 ) ; } } heating { type jouleHeatingSource ; active true ; jouleHeatingSourceCoeffs { anisotropicElectricalConductivity no; sigma 1 .4 2 85 ; } } It can be seen that both source terms have been defined, on the one hand the source term Ser T using the scalaraSemiImplicitSource option type and on the other hand the SJoule Tterm using the jouleHeatingSource option type. The sigma variable represents the electrical conductivity (σ) of the region, which is the inverse of the electrical resistivity previously used. 5.4 Simulation results for BM Once all the above steps are ready, the simulation can be run. To do so, run the command chtMultiRegionSimpleFoam in the terminal. This will take some time to produce the results. The results can be visualised using paraview, as done in the previous steps. In Figures 35 and 36, the results obtained for temperature and voltage are shown respectively. In these figures, it can be seen how this variables evolve across the whole geometry for the case defined. Thermal management study of an electrolytic stack using CFD page. 71 Figure 35: Results obtained for temperature in BM. As can be seen from the results obtained, the simulation tools applied have been appropriate. It can be seen that at the inlet and outlet of the BM, the temperature remains constant at 700ºC. Moreover, the source term due to the electrolysis reaction has been properly introduced, since the region with the highest temperature is the MEA, specifically the centre. This is because the heat generated is extracted from the uppermost layers of the MEA, causing the temperature in the centre of the MEA to rise because the conditions set do not allow for the evacuation of all the internal heat generated. That is, the heat produced by the MEA is extracted through the BPs, so the temperature of the MEA does not increase as much in the outer layers of the BPs. On the other hand, this temperature increase in the MEA is not only due to the heat generated by the electrolysis reaction, but also to the Joule effect. As can be seen in Figure 36, the voltage conditions have been correctly defined. The voltage drops from the MEA, since it is the region that opposes resistance to the passage of electric current. This electrical resistance (which is much higher than that of the BPs) causes heat to be generated by the Joule effect at the MEA. However, this heat produced is very low compared with the heat produced by the operating temperature. That is why the temperature of the BM does not rise so much. page. 72 Report Figure 36: Results obtained for voltage in BM. With this BM a first approach to the electrolyzer has been made, and it has been possible to verify that the simulation techniques that were intended to be used are correct. Therefore, in the following section we will apply a similar proceedure to realize the SOEC model. Thermal management study of an electrolytic stack using CFD page. 79 (sigma) and voltage are defined, as these regions will not have source terms. In addition to the pressure and temperature files (which are also defined for solid regions), it will be necessary to add velocity (U) and hydrostatic pressure (p_rgh) files. An example of the velocity boundary condition file for the cathode channels is shown below. FoamFile { version 2 . 0 ; format a s c i i ; class volVectorField ; obje c t U; } // ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗// dimensions [0 1 −1 0 0 0 0 ] ; in te rn al F i e l d uniform (2.305 0 0) ; boundaryField { inlet_anode { type fixedValue ; value uniform (2.305 0 0) ; } " .∗" { type zeroGradient ; } } Finally, it will be necessary to calculate the value of the source term Ser T, which will be imposed as a generation constant in the MEA. To do this, it is necessary to calculate what is the current flowing through the SOEC when a voltage of 1.3 V is imposed (as in the BM). Following the same procedure as the one used before to obtain the thermophysical properties, the electrical resistivity for the MEA after homogenization is obtained. Table 7: Electrical resistivities for each MEA layer. Layer % Resistivity Hydrogen electrode 81 23.44·10−3 Electrolyte 2.4 23 Intermediate layer 2.4 27.5 Oxygen electrode 14.2 2.6·10−9 page. 80 Report This resistivity obtained is 1.23Ωm. It is obtained by using Table 7and Equation 74. Continuing with the same procedure as above (Section 29), the electrical resistance of the MEA can be calculated with the resistivity obtained and the length and area of the MEA. As seen above in the BM, the electrical resistance of the BPs can be neglected as it is very low compared to that of the MEA. This is why the equivalent resistance of the SOEC will be assumed to be that of the MEA, as shown below: RMEA = 1.23Ωm0.168 ·10−3m 8.1·10−3m2= 0.026Ω Req =  RBP 1+  RBP 2+RMEA =RMEA = 0.026Ω As can be seen, the electrical resistance obtained is much lower than that obtained previously. This is due to the fact that the thickness of the MEA is very low, since this layer for this model is 0.168mm specifically. Finally, applying Ohm’s law, the value of Iis obtained. This value obtained must be expressed in A/m2, so the value obtained is divided by the area through which the current passes (A= 8.1·10−3m2). V=I·Requ →1.3V=I·0.026Ω →I=1.3V 0.026Ω = 50A iz=I A=50 8.1·10−3m2= 6172.84A/m2 With this data and knowing that the value of T·∆Sfor 700°C is 55kJ/mol, the Faraday constant, F, is 96485C/mol and the value of ∆RE is 10.168 ·10−3m, the value of the source term due to the electrolysis reaction can be obtained as shown below. Ser T=izT∆S ∆RE2F=6172.84A/m2·55,000J/mol 96485C/mol ·10.168 ·10−3·2= 173030.39W/m3 This value may seem very high, but it should not be forgotten that it is expressed in cubic meters. Multiplying it by the volume of the MEA region gives a value of 14.25Wfor the source term Ser T. Furthermore, as already mentioned, this value is higher than that obtained for the BM due to the fact that the electrical current flowing is much higher compared to that of the BM, since the electrical resistance of the device is much lower. Thermal management study of an electrolytic stack using CFD page. 81 In the fvOptions file, the Joule effect sources and heat emitted by the electrolysis reaction are defined. An example of the defined file is shown below. In the case of the source term due to the electrolysis reaction, in volumeMode it is defined as absolute, since the value obtained after multiplying by the volume of the region will be defined. Also in this file, the electrical conductivity (the inverse of the electrical resistivity) of the MEA is defined. FoamFile { version 2 . 0 ; format a s c i i ; cl as s dictionary ; location " system " ; obje c t fvOptions ; } // ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗// MEA_AbsoluteEnergySource { type scalarSemiImplicitSource ; active true ; selectionMode cellZone ; cellZone MEA; volumeMode absolute ; sources { h ( 14.25 0 ) ; } } heating { type jouleHeatingSource ; active true ; jouleHeatingSourceCoeffs { anisotropicElectricalConductivity no; sigma 0 . 8 1 3 ; } } 6.4 Simulation results for SOEC Once all the above steps are ready, the simulation can be run. To do so, run the command chtMultiRegionSimpleFoam in the terminal. As there are now two regions that are fluid, the computational cost will be much higher and therefore the simulation will take longer to produce page. 82 Report results. The results can be visualised using paraview, as done in the previous steps. Figure 42 shows the temperature results obtained for the SOEC. Figure 42: Temperature results obtained for SOEC. The results obtained show the expected behaviour, which was already observed in the BM. Figure 42 shows how the part that reaches the highest temperature is the MEA region. This is because it is in this region, in addition to being in contact with the cathode and anode fluids, that the two source terms explained above are present. These source terms, despite being of little relevance with respect to the heat contributed by the fluids, cause the temperature of the MEA to increase a little more compared to the rest of the solid regions. This heat generated in the MEA is evacuated to the outer regions. This is why in the parts of the BPs that are in contact with the MEA, a considerable temperature increase is also observed. This heat generated in the MEA is evacuated to the outer regions. This is why in the parts of the BPs that are in contact with the MEA, a considerable temperature increase is also observed. In Figure 43, the temperature map of two cross sections of the SOEC can be seen. Here it can be seen in more detail that the behaviour described above takes place at both the anode and the cathode. Furthermore, it can be seen in this figure that the highest temperature is reached in the part of the MEA that is in contact with the anode and cathode channels. Thermal management study of an electrolytic stack using CFD page. 83 Figure 43: Temperature map of the SOEC cross section. Figure 44: Voltage results obtained for SOEC. page. 84 Report On the other hand, Figure 44 shows the voltage results obtained for the SOEC. In this case, it can be seen that the BP1 region is at 1.3V while the BP2 region is at 0V. The voltage drop occurs at the MEA, which is where the Joule effect is relevant. This is because the MEA has been modelled as the resistive part, as opposed to the BPs, where the electrical resistance is negligible. Figure 45 shows this potential drop in the MEA in more detail. However, as with the temperature results, the size of the MEA compared to the BPS is very small, which is why the results do not perfectly capture what is happening in this region. Figure 45: Voltage map of the SOEC cross section. 6.5 Future steps As mentioned throughout this document, the scope of this work was to achieve a first approximation of a thermal analysis of a SOEC type electrolyzer by means of CFD. This is why the results obtained are within the expected range. However, it is clear that there are points for improvement that could be carried out in the future in order to achieve a realistic representation of the heat management of an electrolyzer. One of these points would be to eliminate some of the assumptions that have been made in this project. For example, the MEA could be modelled as a porous medium and not as a solid, since in reality some of its layers are porous. Another simplification to be modelled would be to use the correlations obtained for each property of each material, and to calculate their values at each time step. Moreover, it could also be interesting to model the transformation of steam to hydrogen at the cathode, and of air to its increased oxygen content at the anode. That is, to represent the electrolysis reaction itself and the change of species content in the anode and Thermal management study of an electrolytic stack using CFD page. 85 cathode channels. As for the mesh, the refinement of the mesh could be improved in order to capture the heat management in more detail. A scaling of the MEA region could also be done in order to get a more detailed picture of what is happening in that region, since in the results obtained it looks like a very thin region compared to others. In addition, a set of SOECs could also be simulated, which is known as a stack, since these devices are stacked when used on an industrial scale, with up to 100 SOECs in the same stack. Finally, it should not be forgotten that these CFD projects are done with computers that have a high computational capacity. The project, on the other hand, has been developed with a personal computer that is far from these characteristics. In the future, the following parts could be developed with a more powerful computer, which would allow to obtain more accurate results and reduce the time of each simulation. page. 86 Report 7 Economic assessment This section is dedicated to the evaluation of the economic cost of the project. Since the entire project was carried out in one place, the costs have been divided into human resources costs and energy costs. 7.1 Human resources cost Human resources cost is understood as those costs associated to the person/people that has work on creating the present project. It is related with the time spent by the student in the tasks that have been necessary to carry out the project, such as for example, the learning of the necessary knowledge, the development of the fault trees, the drafting of the project report, etc. The economic costs do not include the revision tasks performed by the tutors, so only the cost of the student has been used. The cost of an engineer with less than 5 years of experience has been used as the cost of the student. Table 8: Human resources cost summary Personnel Cost/hour [€/h] Spent time [h] Total cost [€] Student 38 300 11400 Total 11400 (1) 7.2 Energy cost Energy cost is understood as any expense associated with the energy consumption o the equipment used throughout the project. As it is a purely office job, this cost is fundamentally for the electricity consumed by the computer equipment. To calculate this cost, it will be used as data the energy consumption of the equipment, the average price of electricity per kWh in Spain in 2023 and the working hours, as shown in Equation 76. Energy cost = Consumption [kW] ·Working hours [h] ·Electricity cost [€/kWh] (76) The average price per kWh in Spain in 2023 was 0.2064€[26]. In Table 9, the cost associated to each equipment and the total cost for energy is obtained. Table 9: Energy cost summary Equipment Consumption [kW] Working hours [h] Electricity cost [€] Screen 0.019 300 1.17 Computer 0.065 300 4.02 Total 5.19 (2) Thermal management study of an electrolytic stack using CFD page. 87 Once the human resources and energy costs have been calculated, it can be determined which is the total cost of the project, applying the corresponding taxes to both parts. Results are shown in table 10. Table 10: Final costs Concept Calculus applied Total [€] Sum of costs (1) + (2) 11405.19 (3) Surplus for contingencies (15%) (3) ·0.15 1710.78 (4) Total without taxes (3) + (4) 13115.97 (5) VAT (21%) (5) ·0.21 2754.35 (6) Total with taxes (5) + (6) 15870.32 As it is shown in Table 10, the cost associated with the elaboration of this final master’s thesis is fifteen thousand eight hundred and seventy point thirty-two euros. page. 88 Report 8 Environmental impact This section evaluates the environmental impact of the development of the project. As this is an office project, only the environmental impact associated to the production of the energy consumed has been taken into account. To asses this environmental impact analysis, CO2emissions and waste production associated with nuclear energy production have been evaluated. The consumption of the equipment is the one used (73.42 kWh) in Section 7.2. A factor provided by the "Comision Nacional de los Mercados y la Competencia"[26] has been used to calculate CO2 emissions and waste production. Results are shown in Table 11. Table 11: Results for CO2emissions and nuclear waste produced Type Pollutant Factor Total generated Emissions Carbon Dioxide (CO2) 0.25 kg/kWh 18.35 kg Nuclear waste High radioactive waste 0.49 mg/kWh 35.97 mg