Lukas Kasper

Kasper

**HYBRID STORAGE USING THE FINITE ELEMENT METHOD**

**MODELING OF THE PHASE CHANGE MATERIAL OF A** 

## **MODELING OF THE PHASE CHANGE MATERIAL OF A HYBRID STORAGE USING THE FINITE ELEMENT METHOD**

Lukas Kasper MODELING OF THE PHASE CHANGE MATERIAL OF A HYBRID STORAGE USING THE FINITE ELEMENT METHOD

This publication emerged from the project HyStEPs, which is funded by the Austrian Research Promotion Agency (FFG) with grant number 868842. It was supported by the Institute for Energy Systems and Thermodynamics (IET) of the Faculty of Mechanical and Industrial Engineering, TU Wien.

Lukas Kasper

## **MODELING OF THE PHASE CHANGE MATERIAL OF A HYBRID STORAGE USING THE FINITE ELEMENT METHOD**

Cite as:

Kasper, L. (2020). *Modeling of the Phase Change Material of a Hybrid Storage using the Finite Element Method.* TU Wien Academic Press. https://doi.org/10.34727/2020/isbn.978-3-85448-037-2

#### **TU Wien Academic Press 2020**

c/o TU Wien Bibliothek TU Wien Resselgasse 4, 1040 Wien academicpress@tuwien.ac.at www.tuwien.at/academicpress

ISBN (print): 978-3-85448-036-5 ISBN (online): 978-3-85448-037-2

Available online: https://doi.org/10.34727/2020/isbn.978-3-85448-037-2

(CC BY-SA 4.0). https://creativecommons.org/licenses/by-sa/4.0/

Media proprietor: TU Wien, Karlsplatz 13, 1040 Wien Publisher: TU Wien Academic Press Author (responsible for the content): Lukas Kasper Production: Facultas Verlags- und Buchhandels AG

Proofreading: Jennifer Weck

This publication is a revision of the thesis *Modeling of the Phase Change Material of a Hybrid Storage using the Finite Element Method*, written by Lukas Kasper at TU Wien on June 12, 2019. This thesis was published at TU Wien: https://permalink.catalogplus. tuwien.at/AC15369397.

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License

It was written in the framework of the Master programme *Physical Energy and Measurement Engineering* at the Department of Physics of TU Wien. The thesis was supervised and graded by Univ.-Prof. Dr. René Hofmann (supervisor and reviewer), Dr. Alexander Schirrer (co-supervisor) and Dr. Martin Koller.

## Preface

When Lukas Kasper reached out to me and my colleague Priv.-Doz. Dr. techn. Alexander Schirrer in 2018, looking for an interesting and application-oriented topic for his diploma thesis, we seized the opportunity to include him in our recently started research project, Hybrid Storage for Ecient Processes [HyStEPs]. The aim of HyStEPs is to develop an all-in-one storage concept to increase the eciency and flexibility of industrial processes. This project provides interesting research questions and work items that could be answered in the course of a thesis under qualified supervision. The following thesis is the result of the successful symbiosis in the course of this project between an engaged student and experienced scientists. At the TU Wien, great emphasis is put on the participation of students in current research activities. At the Institute for Energy Systems and Thermodynamics, we try to involve students in our research groups as early as possible.

Experience has shown that the concept of research-oriented teaching provides considerable advantages for students, researchers and teachers. Students can benefit from research engagement through increased motivation, a great work environment and direct interchange with colleagues with professional expertise. Students have the unique opportunity to obtain the latest findings from research in teaching and can gain valuable insights into ongoing work. On the one hand, this can lead to a broader range of methodical results; on the other hand, it can also encourage application in a wider variety of subject areas. Methodically, there is the option to deepen theoretical or analytical investigations in related fields. Finally, diploma theses are a viable, compact format for collaboration with industrial partners. These aspects not only enable fruitful research work - they also raise the quality of teaching, with the goal to develop scientific excellence and impart comprehensive competence. This diploma thesis serves as an excellent example in that regard.

Lukas Kasper is currently working on his doctoral thesis, wherein he continues his research in the area of modelling and optimization of industrial energy systems.

Univ.Prof. Dr. techn. Ren´e Hofmann

January, 2020

## Acknowledgements

I am very grateful to my thesis advisor, Prof. Ren´e Hofmann, who gave me the chance to engage myself in this research activity. It was a great experience to be involved in this challenging topic. I wish to thank him for his invaluable advice, which has been essential for this work.

Thank you to all my colleagues from the Institute for Energy Systems and Thermodynamics for the fruitful teamwork and for involving me in the HyStEPs project. I am particularly grateful for the assistance given by Dr. Martin Koller, who also helped me proofread this thesis.

Furthermore, I want to acknowledge my colleagues from the Control and Process Automation Research Unit, especially for conducting numerous simulations on their spontaneously arranged computing grid.

Special thanks to Dr. Alexander Schirrer for lending his expertise in MATLAB and computer science in general.

Finally, I would like to express my gratitude to my parents and my girlfriend Nicki for their constant support throughout my studies. Thank you for your continuous encouragement and for your unfailing patience.

## Abstract

The progressive decarbonisation process and the increased share of renewable energy sources in the grid has increased the need for the development of new approaches to store energy and to increase the eciency of industrial processes. Increasingly, latent heat thermal energy storage systems are being used, which exploit phase change materials to store a large amount of energy during the phase change. A novel approach to increasing the eciency of the commonly used Ruths steam storage is currently being investigated in the project HyStEPs, a project funded by the Austrian Research Promotion Agency (FFG) with grant number 868842. In this concept, a container filled with phase change material is placed at the shell surface of the Ruths steam storage.

In this diploma thesis, and in contribution to the HyStEPs project, the phase change material of this hybrid storage is modelled in two dimensions using the finite element method. The apparent heat capacity method is applied in a MATLAB implementation and considers heat transfer by both conduction and natural convection. Furthermore, the developed code can handle any desired layout of materials arranged on a rectangular domain. The model was successfully validated using an analytical solution and experimental data, and a cross-validation showed corroboration with the results of the CFD software ANSYS Fluent. A parameter study was conducted and the behaviour of di↵erent dimensions and orientations of the phase change material cavity was investigated. The e↵ect of natural convection was found to cause significantly varying behaviour of the studied cavities with di↵erent orientation during the charging process, while the e↵ect was found to be negligible during the discharging process.

## Table of Contents



## Nomenclature

### Variables and Parameters



#### Abbreviations


## List of Figures




## List of Tables


## List of code Listings


# Chapter 1 Introduction

## 1.1 Motivation

Climate change and the management of the earth's resources are two major challenges of our time. The progressive decarbonisation process and the increased share of renewable energy sources in the grid has increased the future demand for thermal energy storage in existing industrial plants and processes is required in the future. The eciency of industrial processes could be increased by balancing steam production and consumption. The Ruths steam accumulator is a well-established thermal storage technology in a variety of industries, including the steel industry.

In order to increase the storage capacity of such steam accumulators, an innovative hybrid storage concept is currently being developed by the HyStEPs project (Hybrid storage for ecient processes), to which this diploma thesis makes a contribution. HyStEPs is a project funded by the Austrian Research Promotion Agency (FFG) with grant number 868842, as part of the Climate and Energy Fund KLI.EN. The project includes the project coordinator Austrian Institute of Technology (AIT) and the partners TU Wien, voestalpine Stahl Donawitz GmbH and EDTMAYER Systemtechnik GmbH.

The basic approach of HyStEPs is to encase Ruths steam storages in innovative latent heat accumulators. Furthermore, a state of charge measurement for the hybrid storage tank should ensure optimal control and operation of both storage types during charging and discharging. This innovative solution is based on an AIT patent application.

### 1.2 Problem statement

Recently, the applications for latent heat thermal energy storages have increased substantially. The main advantage of these thermal energy storages is the high storage density and the ability to store energy at a nearly constant temperature level, using phase change materials (PCMs). This reduces thermal losses compared to sensible heat storages, since energy can be stored at a lower temperature levels. A large amount of energy is absorbed or released during the phase change of such materials, most commonly between the liquid and solid state. The main disadvantage is the low thermal conductivity of most PCMs, which necessitates various methods to increase the heat transfer between the PCM and the heat transfer medium. One approach is the addition of fins inside the PCM containing cavities built from materials with high thermal conductivity, such as aluminium.

In the HyStEPs project, a number of di↵erent design options for encasing cylindric Ruths steam accumulator with a latent heat material have been investigated and are still open for discussion. The most promising of these options, and one that is relatively easy to build, is that of PCM containers with aluminium fins in the plane of the cylinder axis. A schematic illustration of this hybrid storage concept is shown in Figure 1.1. It should be noted that a design optimization of hybrid storage concepts for industrial applications was developed at the same time this thesis was conducted. This optimization was published by fellow HyStEPs project researchers Hofmann et al. [1].

Since heat is transferred in the PCM by both heat conduction and natural convection in the liquid regions, the development of the melting front in the latent heat storage cavity can take very complex forms. This leads to a rather dicult characterisation of the charging and discharging behaviour of the hybrid storage. To determine the position(s) of the melting front(s) in the cavity, a model must be developed which helps to finalize decisions regarding dimensioning and placement of sensors.

*Energies* **2019**, *12*, 898 18 of 25

PCM modules. The implementation of the theoretical ideal design proposal into the direction to reality must be carried out step by step from the calculation to the functional model and further from the test bench to the prototype in compliance with the corresponding safety regulations and must comply with all regulations of the existing strength guidelines, as well as the construction and operating

LHTES modules [43]. Figure 1.1: Schematic sketch of the hybrid storage design consisting of the Ruths steam storage and the surrounding PCM modules [2]

**Figure 6.** Schematic sketch of the hybrid storage design consisting of the RSS and the surrounding

### 1.3 Aim of this work

view A-A. According to [43].

examinations, and inspections.

As mentioned above, one of the challenges in the HyStEPs project is the complex charging and discharging behaviour of the PCM material in the cavities mounted to the Ruths steam storage. As contribution to this project, a numerical model is developed in the scope of this diploma thesis to simulate the evolution of melting and solidification during possible operation scenarios.

The purpose of this model is, on the one hand, the prediction of the exact behaviour of the latent heat storage for di↵erent operation scenarios and design options. This aids the development of the final HyStEPs storage construction, as well as the choice of PCM material, which should also have optimal characteristics in regards to the operational cycle of the steam storage.

On the other hand, the model should provide a basis for developing a reduced order model, which is part of the operation and control strategy of the HyStEPs project and allows full exploitation of the potential of the latent heat storage. The approach for this purpose is a model predictive control based on state of charge measurements in the PCM. Thus, optimal sensor placement must be evaluated, which could be made easier by using a reduced model developed on the basis of the fundamental model created in this thesis.

In this work, a preliminary parameter study was conducted to estimate the charging/discharging speed of di↵erent PCM cavity designs and to study the e↵ect of natural convection in this specific application.

In comparison to commercially available code, the MATLAB model developed should o↵er high usability in regards to running and evaluating parameter studies, which can be run parallel on any desired number of computing units. Furthermore, modifications to the existing code should be easy to make, aiding model reduction by using MATLAB directly.

## 1.4 Methodological approach

The stated problem was analysed with regard to mathematical description, approximation possibilities and numerical modelling options. To this end, a thorough study of relevant academic literature was undertaken. A good review of these topics with regard to phase change materials can be found here [3].

The chosen discretization method of the continuous model is the finite element method, which is well established in applications in structural mechanics, as well as fluid- and thermodynamics. After setting up the fundamental model equations, using temperature as the dependent variable and the apparent heat capacity method to account for the latent heat of phase change, the finite element approximation was derived.

The modelling of the problem stated in Section 1.2 was done in MATLAB [4], which presented a numerical computing environment with very ecient matrix manipulation operations. This is especially suitable for the finite element method.

The model was then validated using an analytical solution for the heat conduction problem and experimental data taken from literature for the problem considering natural convection.

Furthermore, to investigate di↵erent aspects of the latent heat storage cavity of the HyStEPs hybrid storage, simulations were carried out in which certain geometry paramaters are varied and their e↵ect on overall behaviour is analysed. Thus, qualitative statements can be made and the here used methodology used can act as a reference for future investigations.

### 1.5 Structure of this thesis

After this short introduction into the thematics of the problem, this diploma thesis is structured as follows.

In Chapter 2, the underlying theoretical framework for the development of a numerical model describing phase change problems is given. This summary, based on literature review, should provide the reader of this thesis with the most essential knowledge to comprehend the documentation in the subsequent chapters.

The assumptions for the development of the computational model are explained in Chapter 3, where the discretized equations and boundary conditions are also derived. These were later implemented in the form of a MATLAB model, which is documented in Chapter 4. Therein, the basic structure of the code and implementation approaches are explained, and Section 4.3 acts as a guide on how to handle the final version of the code. Detailed instructions on how to specify simulation parameters are available, facilitating use of the code by unfamiliar users.

Chapter 5 deals with the validation of the developed model, which is presented in three steps: The comparison of simulation results considering only heat conduction with the analytical solution of the Stefan problem, the validation of the convection model by experimental data [5], and the cross-validation of the full model with the CFD software ANSYS Fluent. In this chapter of the thesis, very good agreement with the comparative data was obtained and the developed model was therefore successfully validated. In addition, simulations were carried out for di↵erent values of the computational parameters and subsequently analysed. The e↵ect of the chosen mesh size, time step size and mushy region temperature range are discussed, as are possible complications and how to avoid them.

Chapter 6 is a parameter study of the HyStEPs PCM cavity. Within the scope of this thesis, a variation of di↵erent cavity dimensions was simulated and their characteristics evaluated. Furthermore, the e↵ect of natural convection was investigated for di↵erent orientations with respect to position on the steam storage.

The diploma thesis is concluded in Chapter 7, which includes a short summary of the conducted tasks and the obtained results. The scope and limitations of the developed MATLAB model, and of the assertions gained thereby, are discussed in Section 7.2. Finally, improvements and future research objectives are suggested in Section 7.3.

## Chapter 2

## Theoretical framework

### 2.1 Basics of thermodynamics

#### 2.1.1 Principles of thermal energy storage systems

Thermal energy can be stored as a change in the internal energy of a material as either sensible heat, latent heat or thermochemical heat, or a combination of these. Sensible heat storage utilizes the heat capacity and the temperature di↵erence of the material during the process of charging and discharging. The amount of heat stored

$$Q = \int\_{T\_i}^{T\_f} mc\_p dT\tag{2.1}$$

depends on the specific heat capacity *c<sup>p</sup>* at constant pressure *p* of the material, the temperature di↵erence (*T<sup>f</sup> Ti*) and the mass of the storage material *m*, under the assumption of constant *cp*.

Latent heat storage systems are based on the absorption or release of thermal energy when the storage phase change material (PCM) undergoes a phase change. The heat storage capacity is then given by

$$Q = \int\_{T\_i}^{T\_m} mc\_{p,solid} \, dT + ma\_m \Delta l\_m + \int\_{T\_m}^{T\_f} a\_m mc\_{p,liquid} \, dT \,, \tag{2.2}$$

with the specific latent heat *l<sup>m</sup>* and the liquid fraction of the total mass *a<sup>m</sup>* [6]. The specific heat capacity at constant pressure *c<sup>p</sup>* will simply be denoted as *c*, or also as *csolid* for the constant heat capacity of a solid medium and as *cliquid* for the constant heat capacity of a liquid medium.

Phase transitions occur in di↵erent forms:

Solid-solid transitions occur where energy is stored as a material is transformed from one crystalline phase to another. These transitions generally have small latent heat and small volume changes and o↵er the advantages of less stringent container requirements and greater design flexibility [7].


The phase change from solid to liquid (liquefaction) or liquid to solid (solidification) is preferred for the use case in this thesis because the operating pressure is lower than that of liquid to gas or solid to gas phase changes. From here on, this thesis discusses exclusively the solid to liquid or liquid to solid type of phase change.

Liquefaction and solidification can either occur at a constant melting temperature (isothermal phase change) or in a temperature range around the melting temperature *Tm*, which is often called a mushy phase change. Numerous materials, especially pure crystalline substances and eutectics show a sharply defined value for the melting temperature *T<sup>m</sup>* [8], but many PCM materials show a significant melting range. This aids some numerical models in avoiding discontinuities, which can cause problems, as we will see later in this thesis. The change of enthalpy during the phase change for these two types is illustrated in Figure 2.1.

#### 2.1.2 Heat transfer mechanics

In a solid medium, heat is transferred via heat conduction by microscopically colliding particles, including molecules, atoms and electrons, which shift their internal energy from one particle to another. The rate of the energy transfer can be shown as a

Figure 2.1: Enthalpy *H* versus temperature *T* with (a) mushy phase change and (b) isothermal phase change [9]

*Note. L* denotes the specific latent heat.

function of the temperature di↵erence (temperature gradient) in the medium and its conductive properties. This basic law, also known as Fourier's law, is expressed as

$$\mathbf{q} = -k\nabla T\,,\tag{2.3}$$

where ~*q* is the thermal flux and *k* is the thermal conductivity of the medium [10]. In a liquid medium heat is transferred both by conduction and convection. Convection is the transfer of thermal energy via the movement of fluid, and is usually driven by pressure di↵erences in the liquid. In the absence of outer forces, so-called natural convection currents occur, due to temperature-triggered density variations. This leads to a thermal flux

$$\mathbf{q} = \rho c T \mathbf{u} \; , \tag{2.4}$$

where ⇢ is the density of the liquid and u is the velocity of the liquid.

The energy equation for transient heat transfer by both conduction and convection [11] is

$$\frac{\partial(\rho cT)}{\partial t} = \nabla(k\nabla T) - \nabla(\rho cT \cdot \mathbf{u})\ . \tag{2.5}$$

According to the literature, this equation is also called the convection-di↵usion equation or advection-di↵usion equation, depending on the context [12]. Aside from the conservation of energy (2.5), conservation of mass and momentum must be fulfilled.

$$\frac{\partial(\rho)}{\partial t} = \nabla(\rho \mathbf{u})\tag{2.6}$$

For an incompressible fluid with @(⇢) @*<sup>t</sup>* = 0, this continuity equation takes the form

$$
\nabla(\rho \mathbf{u}) = 0 \; . \tag{2.7}
$$

The momentum balance, also known as the Navier-Stokes-equation, takes the following form for an incompressible Newtonian fluid with constant material properties *µ* = *const.,* ⇢ = ⇢<sup>0</sup> = *const.* :

$$
\rho\_0 \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) - \nabla \left( \mu \nabla \mathbf{u} \right) + \nabla p = \mathbf{f} \tag{2.8}
$$

including the pressure *p*, the dynamic viscosity *µ* of the liquid and the external force density f, which usually only consists of the density and the gravitational acceleration g.

#### 2.1.3 Boundary conditions

To solve a system of di↵erential equations, boundary conditions must be established, as must an initial condition for the solution. The main types of boundary conditions used in thermodynamics for the temperature *T*(⌦*, t*) at a boundary ⌦ of the domain ⌦ are listed below [14].

Dirichlet boundary conditions prescribe a temperature at the boundary:

$$T(\delta\Omega, t) = T\_{boundary}(t) \,. \tag{2.9}$$

Neumann boundary conditions prescribe a heat flux q*boundary*(*t*) at the boundary. By using Fourier's law (2.3) these can be written as

$$-k\frac{\partial T}{\partial \mathbf{n}}\left|\_{\delta\Omega,t} = \mathbf{q}\_{boundary}(t)\right|,\tag{2.10}$$

where n is the outward unit normal to the boundary surface [9]. For perfectly insulated boundaries (adiabatic wall) this means

$$-k\frac{\partial T}{\partial \mathbf{n}}\mid\_{\delta\Omega, t} = 0\ . \tag{2.11}$$

Robin boundary conditions prescribe a heat exchange with the surrounding boundary at the ambient temperature *Tboundary*(*t*) with a specified thermal conductivity ↵*boundary*:

$$-k\nabla T \mid\_{\delta\Omega, t} = \alpha\_{boundary} \cdot \left( T \mid\_{\delta\Omega, t} - T\_{boundary}(t) \right) \,. \tag{2.12}$$

This type of boundary condition is often also called Fourier condition or convective boundary condition [9].

It should be noted that there are more types of boundary conditions which are not relevant for this thesis and will therefore not be discussed. One example would be radiation boundary conditions when considering heat radiation e↵ects or periodic boundary conditions, which are useful for symmetric problems. Combinations of the boundary condition types mentioned are sometimes used.

The most common type of boundary condition for the velocity u(⌦*, t*) is the no-slip condition

$$\mathbf{u} \mid\_{\delta\Omega, t} = 0 \; \text{s} \; , \tag{2.13}$$

which essentially prescribes that the fluid at the boundary has zero velocity relative to the boundary, which implies zero velocity in case of a resting boundary. However, for low pressure or high velocity, as well as for low viscosity fluids, a full-slip or partial-slip condition could prove more accurate [13].

#### 2.1.4 The Stefan problem

The analytical solution of the system of di↵erential equations given by (2.5), (2.7) and (2.8) in combination with certain initial and boundary conditions can only be found for a small class of problems. In general, this boundary value problem must be solved numerically, as discussed in Section 2.2.

The simplest form of a mathematical model describing phase transitions is the classic, one-dimensional Stefan problem. This problem describes the liquefaction of a solid medium, initially at the solidification temperature *T*(*x, t* = 0) = *T<sup>m</sup>* with the boundary condition *T*(*x* = 0*, t*) = *T*<sup>0</sup> *> Tm*. To acquire the analytical solution to the classic Stefan problem, meaning the temperature distribution in the liquid phase and the location of the melting front s(t), the heat equation without convection must be solved. A standard form of the problem can be declared as follows [15]:

The liquid region 0 *x<s*(*t*) @(⇢*cT*(*x, t*)) @*<sup>t</sup>* <sup>=</sup> <sup>r</sup>(*k*r*T*(*x, t*)) heat equation, 0 *<x<s*(*t*), *t >* <sup>0</sup> *T*(*x* = 0*, t*) = *T*<sup>0</sup> *> T<sup>m</sup>* Boundary condition, *t >* 0 *T*(*x, t* = 0) = *T<sup>m</sup>* Initial condition

The free boundary *x* = *s*(*t*)

*lm*⇢ @*s* @*<sup>t</sup>* <sup>=</sup> *<sup>k</sup>* @*T* @*x* Stefan condition *s*(*t* = 0) Initial position of melting front *T*(*s*(*t*)*, t*) = 0 Dirichlet condition at the melting front The solid region *s*(*t*) *<x<* 1 *T*(*x, t*) *t*, *x s*(*t*)

A similarity solution [15] of this problem can be found below, dependent on the Stefan number *St<sup>L</sup>* = *<sup>c</sup>*(*T*0*Tm*) *l<sup>m</sup>* and the thermal di↵usivity ↵*<sup>L</sup>* <sup>=</sup> *<sup>k</sup> c*⇢

$$T(x,t) = T\_0 - \frac{T\_0}{\
erf(\lambda\_{St})} \cdot \operatorname{erf}(\frac{x}{2\sqrt{\alpha\_L t}}) \tag{2.14}$$

$$s(t) = 2\lambda\_{St}\sqrt{\alpha\_L t} \tag{2.15}$$

$$
\lambda\_{St} e^{\lambda\_{St}^2} erf(\lambda\_{St}) = \frac{St\_L}{\sqrt{\pi}}\,. \tag{2.16}
$$

Here, *St* is a parameter defined by the transcendental equation (2.16), dependent on the Stefan number *StL*. This solution for the location of the melting front (2.15) will later be used to validate the developed numerical model of the one-dimensional heat equation.

## 2.2 Numerical methods

To solve the system of di↵erential equations given by (2.5), (2.7) and (2.8), numerical methods to accurately describe the phase change must be implemented. Furthermore, the Stefan problem is complicated by the fact that, in general, more than one moving phase boundary can occur in PCM ([16]).

Di↵erent approaches to this problem have been developed, the most common of which are reviewed in the article [17] on which the following section is based.

### 2.2.1 Enthalpy method

One of the most popular approaches is the enthalpy method which can be illustrated by considering a one-dimensional isothermal phase change problem controlled only by heat conduction. The enthalpy function *H* is defined as the integral of the volumetric heat capacity with respect to temperature:

$$H(T) = \Delta H(T) + \rho f(T)a = \int\_{T\_{ref}}^{T} \rho c(T)dT + \rho f(T)a \tag{2.17}$$

The first term of the enthalpy function (2.17) accounts for the sensible heat of the material, while the second term accounts for the latent heat of phase change. To obtain the total enthalpy for a temperature *T*, the liquid fraction

$$a = \begin{cases} 0 & \text{,solid: } T < T\_m \\ & \text{] } 0, 1 \left[ \begin{array}{c} \text{,} \text{,} \text{mushy: } T = T\_m \\ 1 & \text{,} \text{liquid: } T > T\_m \end{array} \right. \end{cases} \tag{2.18}$$

and the function *f*(*T*), describing the amount of latent heat at the temperature *T*, must be given. This leads to an alternate form of the energy equation (2.5)

$$\frac{\partial \Delta H}{\partial t} = \frac{\partial}{\partial x} \left( \alpha\_L \frac{\partial \Delta H}{\partial x} \right) - \rho a \frac{\partial f}{\partial t} \,, \tag{2.19}$$

with the thermal di↵usivity ↵*<sup>L</sup>* = *<sup>k</sup> c*⇢ .

The enthalpy method leads to a simple phase change problem, since interface conditions must not be taken into account [17].

#### 2.2.2 Apparent heat capacity method

In some applications, it is more suitable to use temperature rather than enthalpy as the dependent variable. In these cases, the apparent heat capacity method is often chosen to model the phase change problem. Here, the governing equation is the energy equation (2.5), but the apparent heat capacity of the material *capp* is modelled to account for the latent heat during the phase change:

$$c\_{app} = \begin{cases} c\_{solid} & \text{,solid: } T\_m - \varepsilon > T \\\\ \frac{\Delta l\_m + c\_{solid} \cdot (T\_m + \varepsilon - T) + c\_{liquid} \cdot (-T\_m + \varepsilon + T)}{2\varepsilon} & \text{,msubs: } T\_m - \varepsilon \le T \le T\_m + \varepsilon \\\\ \end{cases}, \begin{aligned} \text{,solid: } T\_m - \varepsilon > T \le T\_m + \varepsilon \\\\ \text{,liquid: } T > T\_m + \varepsilon \end{aligned} \tag{2.20}$$

This method obviously requires the existence of a phase change temperature range, also called melting range, of the width 2✏ around the melting temperature *Tm*. In modelling phase change problems with small values of ✏, this can lead to problems [11]. In addition, the method is generally not applicable in cases with isothermal phase change. However, it is possible to approximate isothermal phase change problems by introducing an appropriately small artificial mushy region of the width 2✏.

### 2.3 Discretization methods

The basic premise of a discretization method is to replace the original continuous problem by a sequence of finite-dimensional problems [18]. In the past, several methods of discretization have been developed to numerically solve problems of partial di↵erential equations, the most common on which will be explained here. It should be noted that there is a distinction between fixed grid and adaptive grid methods, both of which have advantages for specific problems. In this thesis, only fixed grid methods are used and therefore explained, however, most of the following principles can be employed for adaptive grid methods as well.

In the next sections the main discretization methods with respect to space will be described, although most problems are also continuous in time and therefore also need to be discretized appropriately. Time discretization methods will be outlined in Section 2.3.4.

#### 2.3.1 Finite element method

The first step in any finite element (FEM) simulation is to define a grid of the continuous computational domain. A grid consists of a finite number of non-overlapping elements that cover the whole domain [19]. The most common geometries of elements are triangular elements and quadratic elements, the latter of which are used in the model developed in this work. The connection points between the elements, typically located on the corners of the element, are called nodes. Of course, it is also possible to define elements with more nodes, which enhances the accuracy of the approximation, but requires more computational e↵ort.

Based on the defined grid, or mesh, continuous variables, such as the temperature *T*, can be decomposed by the shape functions [*N*]=[*N*1*, N*2*, N*3*...*] and temperatures at the element nodes *{T}* = *{T*1*, T*2*, T*3*...}*:

$$T = [N] \bullet \{T\} \ . \tag{2.21}$$

The general approach of the finite element method is to multiply the partial di↵erential equation (PDE) with given boundary conditions (the strong form) by a test function, or weighting function, and integrate over the whole simulation domain. Applying Green's first integration theorem (partial integration) leads to a variational formulation, also called the weak form [20]. To discretize this still infinite dimensional equation, the approximation of continuous variables, as in (2.21), is applied to the weak form. In the Galerkin procedure, which will also be used in this thesis, the weighting functions are equal to the shape functions, which leads to a system of equations with the same number of equations as unknowns, that is the number of nodes [19]. In the case of a stationary problem, this is a system of algebraic equations, otherwise this becomes a problem of ordinary di↵erential equations.

Although finite di↵erence (FDM) and finite volume schemes (FVM) are well established in computational fluid dynamics (CFD) applications, including the well-known software package ANSYS Fluent, Nedjar [9] states that especially finite element methods (FEM) are able to handle complex coupled thermomechanical problems with various and complex boundary conditions.

#### 2.3.2 Finite di↵erence method

The finite di↵erence method (FDM) can also use the same grid as in 2.3.1. The basic principle of FDM is that the derivatives in the partial di↵erential equation are approximated by linear combinations of function values at the grid points *xi*. For first order derivatives, the di↵erence quotient of a variable *u*(*x*) is by definition [21]

$$\frac{\partial u}{\partial x}(x) = \lim\_{\Delta x \to 0} \frac{u(x + \Delta x) - u(x)}{\Delta x} = \lim\_{\Delta x \to 0} \frac{u(x) - u(x - \Delta x)}{\Delta x} = $$

$$\lim\_{\Delta x \to 0} \frac{u(x + \Delta x) - u(x - \Delta x)}{2\Delta x} \quad (2.22)$$

for an infinitesimal di↵erence *x* between the grid points. This leads to three types of commonly used approximations:

$$\begin{aligned} \left[\frac{\partial u}{\partial x}\right]\_i &\approx \frac{u\_{i+1} - u\_i}{\Delta x} &\text{forward difference} \\ \left[\frac{\partial u}{\partial x}\right]\_i &\approx \frac{u\_i - u\_{i-1}}{\Delta x} &\text{backward difference} \\ \left[\frac{\partial u}{\partial x}\right]\_i &\approx \frac{u\_{i+1} - u\_{i-1}}{2\Delta x} &\text{central difference} \end{aligned} \tag{2.23}$$

The derivatives, which typically arise in two-dimensional convection-di↵usion problems for the variable *u*(*x, y*) and the corresponding grid point indices *i, j* are calculated for the central di↵erence method, as can be seen below:

$$
\left[\frac{\partial u}{\partial x}\right]\_{i,j} \approx \frac{u\_{i+1,j} - u\_{i-1,j}}{2\Delta x} \tag{2.24}
$$

$$
\left[\frac{\partial u}{\partial y}\right]\_{i,j} \approx \frac{u\_{i,j+1} - u\_{i,j-1}}{2\Delta y} \tag{2.25}
$$

$$
\left[\frac{\partial^2 u}{\partial x^2}\right]\_{i,j} \approx \frac{u\_{i+1,j} - 2u\_{i,j} + u\_{i-1,j}}{4\Delta x^2} \tag{2.26}
$$

$$
\left[\frac{\partial^2 u}{\partial x \partial y}\right]\_{i,j} \approx \frac{u\_{i+1,j+1} - u\_{i+1,j-1} - u\_{i-1,j+1} + u\_{i-1,j-1}}{4\Delta x \Delta y} \,. \tag{2.27}
$$

The derivative of a nonlinear term *u*(*x, y*) *v*(*x, y*) can be calculated via the product rule thus

$$\begin{aligned} \left[\frac{\partial uv}{\partial y}\right]\_{i,j} &= u\_{i,j} \left[\frac{\partial v}{\partial y}\right]\_{i,j} + v\_{i,j} \left[\frac{\partial u}{\partial y}\right]\_{i,j} \approx \\ &\quad \frac{u\_{i,j}}{2\Delta y} (v\_{i,j+1} - v\_{i,j-1}) + \frac{v\_{i,j}}{2\Delta x} (u\_{i+1,j} - u\_{i-1,j}) \, . \end{aligned} \tag{2.28}$$

#### 2.3.3 Finite volume method

Although not implemented in the context of this work, it is worth mentioning another important discretization method, the finite volume method (FVM). As in the finite di↵erence and finite element methods, the first step is the discretization of the computational domain, which in this case means splitting it into non-overlapping elements, or finite volumes. The partial di↵erential equations are transformed into algebraic equations by integrating them over each discrete element, in other words, by taking the integral form of conservation equations and splitting the integrals into small intervals. The finite volume method is strictly conservative, because, since two neighbouring cells share a common interface, the total flux entering a volume is identical to that leaving the adjacent volume [22]. This property makes the FVM the preferred method in computational fluid dynamics problems (CFD) [23].

#### 2.3.4 Time discretization

One of the most used methods for time discretization in numerical problems is the so called ✓-method, in which the time derivative of a variable *T*(*t*) is replaced by a simple di↵erential quotient

$$\frac{\partial T}{\partial t} = \frac{T^{n+1} - T^n}{\Delta t} \,. \tag{2.29}$$

Here, *T <sup>n</sup>* = *T*(*x, tn*) is the variable's value at *t* = *t<sup>n</sup>* and *t* is the time increment according to *tn*+1 = *t<sup>n</sup>* + *t*. The variable *T*(*x, tn*) is usually assumed to be known and used as an initial condition to advance the solution to the next time level. The solution *T* can then be written as

$$T = \theta T^{n+1} + (1 - \theta)T^n \ , \ t\_n \le t \le t\_{n+1} \ . \tag{2.30}$$

The relaxation parameter ✓ 2 [0*,* 1] is most commonly one of the following:

✓ = 1 implicit backwards Euler scheme ✓ = <sup>1</sup> <sup>2</sup> second-order, centered implicit method (Crank-Nicolson-Galerkin) ✓ = 0 explicit forward Euler scheme [19]

It is known that for ✓ *<* <sup>1</sup> <sup>2</sup> , only conditional stability is attained [24]. In the model developed in this thesis, the implicit backwards Euler scheme was chosen.

In Chapter 3.2 the semi-discrete Galerkin-formulation

$$\begin{Bmatrix} [C] \left\{ \dot{T} \right\} + [K] \left\{ T \right\} = [R] \end{Bmatrix} \tag{3.23}$$

(2.31)

will be developed (discrete in space, continuous in time). In this equation, the square brackets denote a matrix and the braces denote a vector of node point values.

If we consider this system of ordinary di↵erential equations and apply the backwards Euler scheme, we evaluate the equation at the time *t* = *tn*+1:

$$\left[ \left[ C \right] \left\{ \dot{T} \right\} \right]\_{n+1} + \left[ K \right] \left\{ T \right\}\_{n+1} = \left[ R \right] \ . \tag{2.32}$$

The time derivative of the temperature vector *{T}* is approximated through the backwards di↵erential quotient

$$\left\{\frac{\partial T}{\partial t}\right\}\_{n+1} = \left\{\dot{T}\right\}\_{n+1} = \frac{\{T\}\_{n+1} - \{T\}\_n}{\Delta t} \,, \tag{2.33}$$

which gives us the value

$$\{T\}\_{n+1} = \{T\}\_n + \Delta t \left\{\dot{T}\right\}\_{n+1} \,. \tag{2.34}$$

Combining equation (2.32) with (2.34) leads, with the expression

$$[M] = [C] + \Delta t \bullet [K] \; , \tag{2.35}$$

to

$$\left[ \dot{M} \right] \left\{ \dot{\hat{T}} \right\}\_{n+1} = \left( \left[ C \right] + \Delta t \bullet \left[ K \right] \right) \left\{ \dot{\hat{T}} \right\}\_{n+1} = \left[ R \right] - \left[ K \right] \left\{ T \right\}\_{n} \ . \tag{2.36}$$

We can therefore derive the value for <sup>n</sup> *T*˙ o *n*+1 as

$$\left\{ \dot{T} \right\}\_{n+1} = \left[M\right]^{-1} \cdot \left( \left[R\right] - \left[K\right] \left\{ T\right\}\_{n} \right) \; , \tag{2.37}$$

if the matrix [*M*] (equation (2.35)) is invertible and nonsingular [20]. The solution for *{T}n*+1 can thus be found through equation (2.34) to be

$$\begin{aligned} \{T\}\_{n+1} &= \{T\}\_n + \Delta t \left\{ \dot{T} \right\}\_{n+1} \\\\ &\{T\}\_n + \Delta t \left[M\right]^{-1} \cdot \left( \left[R\right] - \left[K\right] \left\{T\right\}\_n \right) = \\\\ &\{T\}\_n + \Delta t \left( \left[C\right] + \Delta t \cdot \left[K\right] \right)^{-1} \cdot \left( \left[R\right] - \left[K\right] \left\{T\right\}\_n \right) \ . \end{aligned} \tag{2.38}$$

# Chapter 3 Numerical modelling

In this chapter, the underlying assumptions for developing a computational model of the PCM storage are explained, and the discretized equations and boundary conditions are derived.

Although a number of di↵erent design options are currently being discussed for the PCM storage containers attached to the cylindric Ruth steam storage, one of the most promising and easiest to build calls for PCM containers with aluminium fins in the plane of the cylinder axis. For general fin arrangements, a three-dimensional cavity simulation might be required. However, the geometry of the PCM cavities regarded in this work, in which the depth of the enclosure parallel to the cylinder axis is assumed large enough for wall boundary layer e↵ects to be negligible, suggests symmetry reduction of the simulation problem to single aluminium fin segments and a two-dimensional approach. In the hybrid storage concept in [25] this geometry relates to an aluminium fin structure orientated in the radial-axial plane of an RSS.

This sort of design is modelled in this thesis, although the developed model can be easily modified to suit similar problems. Unfortunately, it is not possible to simulate a design with fins oriented in a plane perpendicular to the cylinder axis. Heat conduction and convection e↵ects would make this a three-dimensional problem and, as such, this would require further extensions to a three-dimensional finite element model.

### 3.1 Model equations and assumptions

#### 3.1.1 Geometry

Although the fins are not exactly parallel, the diameter of the cylindric Ruths steam storage is assumed much larger than the thickness of the attached PCM modules. Therefore, it is reasonable to see the PCM container as rectangular. Also, it is possible that the containers will later be constructed in a rectangular design, due to their simplified construction. The underlying equations of the problem will therefore be developed for a rectangular mesh of finite elements. PCM material, aluminium fins, as well as containment or other materials can later be defined separately for each element.

#### 3.1.2 Material properties

Currently, several di↵erent PCM materials are being discussed for inclusion in the application described in this work. A comparison of di↵erent PCM materials can be found in [26]. Thus far, the most promising material considering the container geometry and operating temperatures of this hybrid storage system is a eutectic mixture of potassium nitrate and sodium nitrate (KNO3-NaNO3). The properties of this mixture, as well as those of a typical containment steel, carbon steel 1.0425, taken from reference [27], are compared in Table 3.1. In this reference, it is noted that despite the theoretical melting temperature of the mixture at 222 C, the onset of melting of the technical grade material used was measured at 219*.*15 C. For the purpose of basic evaluations in this thesis, an approximation of *T<sup>m</sup>* = 220 C is used and is therefore given as such in Table 3.1.

There is no exact number to declare the melting temperature range of (KNO3- NaNO3), but there are results of di↵erent measurements found in literature. In reference [28], the reported result for the onset of melting was 220*.*6 *±* 0*.*3 C and for the onset of crystallization 223*±*0*.*34 C, which would result in a melting range of around 2*.*5 C. It was also noted that the results are in very good agreement with those found in other literature and that neither subcooling nor thermochemical instability have been observed [28].



For most simulations aiming to assess the behaviour of the potential PCM container, a mushy region of 2 C was introduced, otherwise this will be declared. The e↵ect of the mushy region with respect to possible numerical diculties will be discussed later in this thesis.

Also given in Table 3.1 are the well known properties of aluminium and air [29], which are used for certain simulations.

Regarding the heat transfer from the RSS to the PCM containers, the choice of inner heat transfer coecient is not straightforward, since it is highly dependent on the inner thermodynamic state in the RSS and therefore on the operating scenario. If the RSS is charged, condensing steam dominates the heat transfer to the outer wall in the vapour phase, while in the liquid phase, the heat transfer is dominated by liquid water. During discharging, heat transfer through boiling water occurs at the outer wall in the liquid phase, and heat transfer through dry steam occurs in the vapour phase [30]. The respective heat transfer coecients are listed in Table 3.2.

In the parameter studies following in Chapter 6, heat transfer coecients were assumed that are within the range given here. However, the impact of di↵erent heat transfer coecients on the hybrid storage behaviour was not within the scope of this thesis. A further investigation of this research topic was pursued in [31].

Table 3.2: Heat transfer coecients ↵ occurring between the inside of the RSS and its outer wall [32]


#### 3.1.3 Governing equations

Based on the theoretical framework explained in Section 2.1.2, the governing equation for the conduction-convection-model can again be stated as follows:

$$\text{energy equation:} \quad \frac{\partial(\rho cT)}{\partial t} = \nabla(k\nabla T) - \nabla(\rho cT \cdot \mathbf{u}) \tag{2.5}$$

$$\text{continuity equation:} \quad \nabla(\rho \mathbf{u}) = 0 \tag{2.7}$$

$$\text{Navier-Stokes equation:} \quad \rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) - \nabla \left( \mu \nabla \mathbf{u} \right) + \nabla p = \mathbf{f} \tag{2.8}$$

By assuming a constant density ⇢ over the whole operating temperature range and neglecting the time and space dependency of the heat capactiy *c*, which is often found in literature for PCM studies using the apparent heat capacity method [11], the energy equation can be written as

$$
\rho c \frac{\partial (T)}{\partial t} = \nabla (-\mathbf{q}) - \rho c \nabla (T \cdot \mathbf{u}) \ , \tag{3.1}
$$

with the flux q = *k*r*T*. Assuming constant density also means freedom of divergence of the velocity field, thus equation (2.7) is reduced to

$$
\nabla \mathbf{u} = 0 \; . \tag{3.2}
$$

To account for phase change in the material, the apparent heat capacity method is used, as introduced in Section 2.2.2:

$$c\_{app} = \begin{cases} c\_{solid} & \text{,solid: } T\_m - \varepsilon > T \\\\ \frac{\Delta l\_m + c\_{solid} \cdot (T\_m + \varepsilon - T) + c\_{liquid} \cdot (-T\_m + \varepsilon + T)}{2\varepsilon} & \text{,msubs: } T\_m - \varepsilon \le T \le T\_m + \varepsilon \\\\ c\_{liquid} & \text{,liquid: } T > T\_m + \varepsilon \end{cases} \tag{2.20}$$

For better readability, the apparent heat capacity *capp* is henceforth denoted only as *c*.

The Navier-Stokes equation (2.8) cited in the previous chapter is already in a particular form, that of an incompressible Newtonian fluid with constant material properties ⇢*, µ* = *const*. A further approximation is to assume constant pressure and only account for density di↵erences in terms of the force density f. This is known as the Boussinesq approximation, where the buoyancy force

$$\mathbf{f} = -\rho \mathbf{g} = -\rho\_0 \left(\beta \left(T - T\_{ref}\right)\right) \mathbf{g} \tag{3.3}$$

can be calculated by a linear thermal expansion coecient , the reference density ⇢<sup>0</sup> and a reference temperature *Tref* [33]. Thus, the Navier-Stokes equation to be solved is

$$
\rho\_0 \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) - \mu \nabla \left( \nabla \cdot \mathbf{u} \right) + \nabla p = \mathbf{f} \ . \tag{3.4}
$$

#### 3.1.4 Boundary conditions

To successfully solve the system of partial di↵erential equations, certain initial and boundary conditions must be given. An initial condition for the temperature in the simulation domain must be specified, and the values can be fixed for each individual node, if desired. Concerning the domain boundary conditions, the following types have been implemented in the model:

Robin boundary conditions (see (2.12)) have been applied to the east and west domain boundaries:

$$|\mathbf{q}(x=0,y)| = |\mathbf{q}\_{in}(y)| = \alpha\_{in} \bullet (T(x) - T\_{in}) \text{ and} \tag{3.5}$$

$$|\mathbf{q}(x = L\_x, y)| = |\mathbf{q}\_{out}(y)| = \alpha\_{out} \cdot (T(x) - T\_{out})\,\,. \tag{3.6}$$

Adiabatic Neumann boundary conditions have been applied to the south and north domain boundaries:

$$\mathbf{q}(x, y=0) = 0 \text{ and } \tag{3.7}$$

$$\mathbf{q}(x, y = L\_y) = 0 \; . \tag{3.8}$$

The selectable thermal conductivities ↵*in* and ↵*out* and ambient temperatures *Tin* and *Tout* can also be chosen to be dependent of time and x or y respectively. For the sake of better readability, this is not explicitly depicted here. The values of *L<sup>x</sup>* and *L<sup>y</sup>* denote the length of the domain in the corresponding dimension.

Regarding the velocity field u = [*u, v*] *<sup>T</sup>* , no slip boundary conditions are set for the domain boundaries

$$\mathbf{u} = \begin{bmatrix} u, v \end{bmatrix}^T = 0 \ \text{, for } x = L\_x, y = L\_y, x = 0, y = 0 \ \text{,} \tag{3.9}$$

and the velocity is also set to zero if the PCM is solid, as well as for every non-melting material

$$\mathbf{u} = [u, v]^T = 0 \text{ , for } T < T\_m - \epsilon \text{ .} \tag{3.10}$$

### 3.2 Discretization

In the following section, the energy equation (3.1) is discretized using the standard Galerkin finite element method, as introduced in Section 2.3.1. The time discretization implemented later is chosen as the implicit backwards Euler scheme, as explained in section 2.3.4.

The Navier-Stokes equation (3.4) is treated separately via a finite di↵erence discretization, which will be explained in Section 3.3.

#### 3.2.1 Preliminary remarks

The two-dimensional domain (*x* 2 [0*, Lx*], *y* 2 [0*, Ly*]) is split into an equidistant mesh of *nnel* = *n<sup>x</sup> n<sup>y</sup>* finite elements. In this thesis, four-noded bilinear rectangular elements are used, which leads to the total sum of *nnon* = (*n<sup>x</sup>* + 1) (*n<sup>y</sup>* + 1). The dependent variables, meaning the temperature *T* and the velocity u = [*u, v*] *<sup>T</sup>* , can therefore be approximated by their nodal values as

$$T = [N] \cdot \{T\} \; , \tag{3.11}$$

$$u = [N] \cdot \{u\} \text{ and } \tag{3.12}$$

$$v = [N] \cdot \{v\} \; , \tag{3.13}$$

where [*N*]=[*N*1*, N*2*, N*3*, ..., Nnon*] is the vector of shape functions and the braces denote a vector of node point values. In the notation used here, general matrices are represented by square brackets.

The flux can be discretized as

$$\mathbf{q} = -k \cdot [B] \bullet \{T\} \tag{3.14}$$

when using the gradient of the shape function matrix

$$[B] = \left[ \begin{array}{cccc} \frac{\partial N1}{\partial x} & \frac{\partial N2}{\partial x} & \frac{\partial N3}{\partial x} & \dots \\ \frac{\partial N1}{\partial y} & \frac{\partial N2}{\partial y} & \frac{\partial N3}{\partial y} & \dots \\ \end{array} \right] . \tag{3.15}$$

To illustrate the approximation via the four-noded bilinear elements, we consider a finite element in the natural coordinate system ⇠, ⌘ with the length *dx* = 2 and height *dy* = 2 (see Figure 3.1). A function (⇠*,* ⌘) can then be approximated in terms of the nodal values and shape functions as

$$\Phi(\xi,\eta) = N\_1\Phi\_1 + N\_2\Phi\_2 + N\_3\Phi\_3 + N\_4\Phi\_4 \ , \tag{3.16}$$

with the shape functions

$$\begin{aligned} N\_1 &= \frac{1}{4}(1-\xi)(1-\eta) \ , \\ N\_2 &= \frac{1}{4}(1+\xi)(1-\eta) \ , \\ N\_3 &= \frac{1}{4}(1+\xi)(1+\eta) \ , \\ N\_4 &= \frac{1}{4}(1-\xi)(1+\eta) \ [19] \ . \end{aligned} \tag{3.17}$$

Figure 3.1: Four-noded bilinear element in natural coordinate system ⇠*,* ⌘ [19] with node coordinates in brackets

#### 3.2.2 Galerkin's method

Using Galerkin's method on the energy equation (3.1) leads to

$$\int\_{V} \left[ N \right]^T \rho c \frac{\partial T}{\partial t} \, dV + \int\_{V} \left[ N \right]^T \nabla (\mathbf{q}) \, dV + \int\_{V} \left[ N \right]^T \rho c \nabla (T \mathbf{u}) \, dV = 0 \, . \tag{3.18}$$

When applying the divergence theorem R *<sup>V</sup>* <sup>r</sup>*F dV* <sup>=</sup> <sup>R</sup> *<sup>S</sup> F* n *dS* [34] on the third integral we obtain the following relation:

$$\begin{split} \int\_{V} \left[ N \right]^{T} \rho c \nabla (T \cdot u) \, dV + \int\_{V} \nabla ([N]^{T}) \rho c T \mathbf{u} \, dV &= \int\_{V} \rho c \nabla ([N]^{T} \, T \mathbf{u}) \, dV = \\ &\quad \int\_{S} \rho c ([N]^{T} \, T \mathbf{u})^{T} \cdot \mathbf{n} \, dS \, . \end{split} \tag{3.19}$$

Because of boundary condition (3.9), the surface integral over the domain boundary is equal to zero and the third integral of equation (3.18) can be written as:

$$\int\_{V} \left[ N \right]^{T} \rho c \nabla \left( T \mathbf{u} \right) dV = \int\_{V} \nabla \left( \left[ N \right]^{T} \right) \rho c T \mathbf{u} \, dV \, . \tag{3.20}$$

Similarly, inserting the flux q = *k*r*T* into the second integral and applying the divergence theorem leads to

$$\int\_{V} \left[ N \right]^T \nabla \left( \mathbf{q} \right) dV = \int\_{V} \nabla \left[ \left[ N \right]^T \right] k \nabla T \, dV + \int\_{S} \left[ N \right]^T \mathbf{q} \cdot \mathbf{n} \, dS \, \,, \tag{3.21}$$

which is written as

$$\begin{aligned} \int\_{V} \left[ B \right]^{T} k \left[ B \right] \left\{ T \right\} \, dV + \int\_{S1} \left[ N \right]^{T} \alpha\_{out} (\left[ N \right] \left\{ T \right\} - T\_{out}) \, dS + \\ \int\_{S3} \left[ N \right]^{T} \alpha\_{in} (T\_{in} - \left[ N \right] \left\{ T \right\}) \, dS + \int\_{S2} 0 \, dS + \int\_{S4} 0 \, dS \end{aligned} \tag{3.22}$$

when applying the boundary conditions for the boundary *S*, where the index *S*1 denotes the east boundary, *S*3 the west boundary and indices *S*2, *S*4 the south and north boundaries, where zero flux is defined.

Finally, inserting the finite element approximations into equation (3.18) leads to the system of ordinary di↵erential equations

$$\left[ \left[ C \right] \left\{ \dot{T} \right\} \right] + \left[ K \right] \left\{ T \right\} = \left[ R \right] \; , \tag{3.23}$$

with the following system matrices:

$$\left[C\right] = \int\_{V} \left[N\right]^T \rho c \left[N\right] \,dV \,\,, \tag{3.24}$$

$$\begin{aligned} [K] = \int\_{V} \left[ B \right]^{T} \mathbf{k} \left[ B \right] \, dV - \int\_{V} (\left[ B \right]^{T}) \rho c \mathbf{u} \, [N] \, dV \\ \qquad + \int\_{S1} \alpha\_{out} (\left[ N \right]^{T} \left[ N \right] \, dS - \int\_{S3} \alpha\_{in} (\left[ N \right]^{T} \left[ N \right] \, dS) \, , \end{aligned} (3.25)$$

$$\left[R\right] = \int\_{S1} \alpha\_{out} T\_{out} \left[N\right]^T \, dS - \int\_{S3} \alpha\_{in} T\_{in} \left[N\right]^T \, dS \, . \tag{3.26}$$

The system matrix [*R*] is of the dimension *non*⇥1, while the matrices [*C*] and [*K*] are of the dimension *non*⇥*non*, when the velocity u = [*u, v*] *<sup>T</sup>* and the thermal conductivity k = [*kx, ky*] *<sup>T</sup>* = *k* [1*,* 1]*<sup>T</sup>* are evaluated. The time discretization of equation (3.23) is done by application of the implicit backwards Euler method as derived in Section 2.3.4 of this thesis:

$$\begin{aligned} \{T\}\_{n+1} &= \{T\}\_n + \Delta t \left\{ \dot{T} \right\}\_{n+1} \\\\ &\{T\}\_n + \Delta t \left[M\right]^{-1} \bullet \left( \left[R\right] - \left[K\right] \left\{T\right\}\_n \right) = \\\\ &\{T\}\_n + \Delta t \left( \left[C\right] + \Delta t \bullet \left[K\right] \right)^{-1} \bullet \left( \left[R\right] - \left[K\right] \left\{T\right\}\_n \right) \ . \end{aligned} \tag{2.38}$$

#### 3.2.3 Integration scheme

To numerically compute the system matrices derived in the previous section, the integration is split into the finite domains of the defined elements. The integrals are first transformed into the natural coordinate system 1 ⇠ ⌘ 1, which leads to simple integration limits:

$$\int\_{V} f(x, y) \, dV = \int\_{-1}^{1} \int\_{-1}^{1} f(\xi, \eta) \, |\mathbf{J}| \, d\xi d\eta \tag{3.27}$$

The factor *|*J*|* is the determinant of the Jacobian matrix of the transformation J [19], which is defined by

$$
\begin{bmatrix}
\frac{\partial N}{\partial \xi} \\
\frac{\partial N}{\partial \eta} \\
\frac{\partial N}{\partial \eta}
\end{bmatrix} = \begin{bmatrix}
\frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\
\frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta}
\end{bmatrix} \cdot \begin{bmatrix}
\frac{\partial N}{\partial x} \\
\frac{\partial N}{\partial y}
\end{bmatrix} = \mathbf{J} \cdot \begin{bmatrix}
\frac{\partial N}{\partial x} \\
\frac{\partial N}{\partial y}
\end{bmatrix} \cdot \tag{3.28}
$$

Now the numerical integration procedure known as Gaussian quadrature is applied, which approximates a definite integral with a weighted sum over a finite set of points [19]:

$$\int\_{-1}^{1} \int\_{-1}^{1} f(\xi, \eta) \left| \mathbf{J} \right| d\xi d\eta \cong \sum\_{i=1}^{n\_{gp}} \sum\_{j=1}^{n\_{gp}} W\_i W\_j f(\xi\_i, \eta\_j) \left| \mathbf{J} \right|\ . \tag{3.29}$$

*W<sup>i</sup>* are the weight factors, ⇠*<sup>i</sup>* and ⌘*<sup>i</sup>* are the Gauss points and *ngp* is the number of Gauss integration points, which specify that a polynomial of degree 2*ngp* 1 is integrated exactly by the method of Gauss-Legendre [19].

This method, which is the most common form of Gaussian quadrature, is chosen for the integration of the system matrices [*C*], [*K*] and [*R*], with the exception of the velocity containing term

$$-\int\_{V} (\left[B\right]^{T}) \rho c \mathbf{u} \left[N\right] \,dV\tag{3.30}$$

in the system matrix [*K*], where a three-point Gauss-Lobatto quadrature rule is applied. In this rule, two of the three Gauss points lie on the ends of the integration interval. This was found to be more accurate due to the natural inclusion of boundary conditions. This was not further studied in a mathematically rigorous way, but rather, observed through testing. Because of the greater computational expense, this was programmed in such a way, that the term (3.30) is only evaluated for u 6= 0. Gauss points and weight factors of the mentioned integration rules are given in Table 3.3.

Table 3.3: Gauss points and weight factors of the used quadrature rules [35]


### 3.3 Convection modelling

As previously mentioned, the Navier-Stokes equation (3.4) is treated separately from the heat transfer part of the model via a finite di↵erence discretization. In the early stages of this diploma work, an attempt was made to also discretize the Navier-Stokes equation with the finite element method, similar to the published work [11]. However, this method, using a penalty formulation for the pressure, could not be successfully implemented in the MATLAB code, as it did not lead to convergent solutions. To avoid this problem, an open source code, available on the course homepage (http://math.mit.edu/˜gs/cse/) of 'Computational Science and Engineering' at the Massachusetts Institute of Technology, was implemented and extended. The documentation of this code can be found in [36]. In the following sections, the basic functioning of the existing MIT code will be described and later on, the modifications that were made in order to fit the here described problem will be explained.

At this point, it is necessary to point out that the methods in the following sections are merely a summary of the mentioned documentation [36] and not the author's own work, unless declared otherwise.

#### 3.3.1 Two-dimensional Navier-Stokes equations

In two space dimensions, we derive for the nonlinear convective acceleration

$$(\mathbf{u} \cdot \nabla)\mathbf{u} = \left[ \begin{pmatrix} u \\ v \end{pmatrix} \cdot \begin{pmatrix} \partial / \partial x \\ \partial / \partial y \end{pmatrix} \right] \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} \frac{\partial u^2}{\partial x} + \frac{\partial uv}{\partial y} \\ \frac{\partial uv}{\partial x} + \frac{\partial v^2}{\partial y} \end{pmatrix} \tag{3.31}$$

and the term r (r u) can be written as

$$\nabla \left( \nabla \cdot \mathbf{u} \right) = \Delta \mathbf{u} = \begin{pmatrix} \frac{\partial u}{\partial x^2} + \frac{\partial u}{\partial y^2} \\ \frac{\partial v}{\partial x^2} + \frac{\partial v}{\partial y^2} \end{pmatrix} \tag{3.32}$$

The Navier-Stokes equation (3.4) and continuity equation in two space dimensions thus take the form

$$\frac{\partial u}{\partial t} + \frac{1}{\rho\_0} \frac{\partial p}{\partial x} = -\frac{\partial u^2}{\partial x} - \frac{\partial uv}{\partial y} + \frac{\mu}{\rho\_0} \left( \frac{\partial u}{\partial x^2} + \frac{\partial u}{\partial y^2} \right) + \frac{f\_x}{\rho\_0} \text{ and} \tag{3.33}$$

$$\frac{\partial v}{\partial t} + \frac{1}{\rho\_0} \frac{\partial p}{\partial y} = -\frac{\partial v^2}{\partial y} - \frac{\partial uv}{\partial x} + \frac{\mu}{\rho\_0} \left( \frac{\partial v}{\partial x^2} + \frac{\partial v}{\partial y^2} \right) + \frac{f\_y}{\rho\_0} \,, \tag{3.34}$$

$$
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = -\frac{\partial \rho}{\partial t} = 0 \,, \tag{3.35}
$$

whereby *f<sup>x</sup>* and *f<sup>y</sup>* are the components of the buoyancy force (3.3).

#### 3.3.2 Numerical solution approach

The numerical approach to solving the time dependent problem is the fractional step method [37], which splits the Navier-Stokes equations into equations that are significantly simpler to work with. The solution is then found by following a three step approach, which is outlined in the following, where u*<sup>t</sup>* denotes the solution at the time *t* and u*t*+1 the final solution at the time *t* + 1, while u⇤ and u⇤⇤ denote intermediate solutions of the respective split equations.

#### 1. Explicit time step for nonlinear terms

The nonlinear convective acceleration term (3.31) is treated explicitly, which circumvents the implicit solution of a large nonlinear system of equations.

$$u^\* - u\_t - \Delta t \left(\frac{\partial u\_t^2}{\partial x} + \frac{\partial u\_t v\_t}{\partial y}\right) = 0\tag{3.36}$$

$$v^\* - v\_t - \Delta t \left(\frac{\partial u\_t v\_t}{\partial x} + \frac{\partial u\_t^2}{\partial y}\right) = 0\tag{3.37}$$

This, however, introduces a Courant-Friedrichs-Lewy (CFL) condition, which limits the time step by a constant times the spatial resolution [38]. For a firstorder upwind scheme, this CFL condition requires that, for stability, the conditions

$$CFL\_{(x)} \equiv \frac{|u|\_{max} \,\Delta t}{\Delta x} \le 1\tag{3.38}$$

and

$$CFL\_{(y)} \equiv \frac{|v|\_{max} \,\Delta t}{\Delta y} \le 1\tag{3.39}$$

hold [39].

#### 2. Implicit time step for viscosity terms

If the viscosity terms (3.32) were treated explicitly, this would result in a time step restriction proportional to the spatial discretization squared. This part of the equation is therefore handled by an implicit step, at the expense of two linear systems of equations to be solved in each time step:

$$u^{\*\*} - u^\* = \Delta t \frac{\mu}{\rho\_0} \left( \frac{\partial u\_{\*\*}}{\partial x^2} + \frac{\partial u\_{\*\*}}{\partial y^2} \right) + \Delta t \frac{f\_x}{\rho\_0} \text{ and} \tag{3.40}$$

$$v^{\*\*} - v^\* = \Delta t \frac{\mu}{\rho\_0} \left( \frac{\partial v\_{\*\*}}{\partial x^2} + \frac{\partial v\_{\*\*}}{\partial y^2} \right) + \Delta t \frac{f\_y}{\rho\_0} \,. \tag{3.41}$$

#### 3. Pressure correction step

The intermediate velocity field (*u*⇤⇤*, v*⇤⇤) is corrected by the gradient of pressure *pt*+1 to ensure incompressibility.

$$u\_{t+1} - u^{\*\*} = -\frac{\Delta t}{\rho\_0} \frac{\partial p\_{t+1}}{\partial x} \tag{3.42}$$

$$v\_{t+1} - v^{\*\*} = -\frac{\Delta t}{\rho\_0} \frac{\partial p\_{t+1}}{\partial y} \tag{3.43}$$

33

The pressure *pt*+1 is only given implicitly and can be obtained by solving a linear system of equations. Applying the divergence operator r on both sides of equations (3.42) and (3.43), and taking into account the continuity equation ru*t*+1 = 0 leads to:

$$
\Delta p\_{t+1} = \frac{\rho\_0}{\Delta t} \nabla \mathbf{u}^{\*\*} \tag{3.44}
$$

The pressure correction step can therefore be implemented as follows:


When solving the Poisson equation, the question arises on the appropriate choice of boundary conditions for the pressure *p*. In the documentation [36], it is stated that a standard approach is to prescribe homogenous Neumann boundary conditions wherever no-slip boundary conditions are prescribed to the velocity field, which, in this case is on every boundary. This implies that the pressure is only defined up to a constant. However, absolute pressure information is not needed because only the gradient of *p* is needed in the pressure correction step.

It can be noted that combining equations (3.36), (3.40) and (3.42) again yields the full Navier-Stokes equation system (3.33). The same is of course true for the velocity component *v* in equation (3.34).

#### 3.3.3 Spatial discretization

The spatial discretization of the Navier-Stokes equation is done via the finite di↵erence method, as introduced in Section 2.3.2. A staggered grid (see Figure 3.2) based on the finite element grid is used to solve the heat transfer problem, where the pressure *p* is defined in the element centers, *u* is placed in the middle of the vertical element borders and *v* is placed in the middle of the horizontal element borders.

Figure 3.2: Two-dimensional domain discretized into a staggered grid [36]

The second and first derivatives are approximated by centered finite di↵erence stencils, as illustrated here for *ui,j* and ✓@*<sup>u</sup>* @*x* ◆ *i,j* at the grid point (*i, j*):

$$
\Delta u\_{i,j} = \left(\frac{\partial^2 u}{\partial x^2}\right)\_{i,j} + \left(\frac{\partial^2 u}{\partial y^2}\right)\_{i,j} \approx \frac{u\_{i-1,j} - 2u\_{i,j} + u\_{i+1,j}}{\Delta\_x^2} + \frac{u\_{i,j-1} - 2u\_{i,j} + u\_{i,j+1}}{\Delta\_y^2} \tag{3.45}
$$

$$
\left(\frac{\partial u}{\partial x}\right)\_{i,j} \approx \frac{u\_{i+1,j} - u\_{i-1,j}}{2\Delta x} \,. \tag{3.46}
$$

This approximation for the first derivative can yield instabilities [36], but the centered stencil

$$
\left(\frac{\partial u}{\partial x}\right)\_{i+\frac{1}{2},j} \approx \frac{u\_{i+1,j} - u\_{i,j}}{\Delta x} \tag{3.47}
$$

is a stable approximation for @*<sup>u</sup>* @*<sup>x</sup>* between two grid points [36]. This aids the staggered grid method, where this position is exactly the position of the pressure *pi,j* .

For the nonlinear terms, a combination of the centered di↵erencing and upwinding approach is implemented, with a transition parameter 2 [0*,* 1] defined as

$$\gamma = \min\left(1.2 \cdot \Delta t \cdot \max\left(\max\_{i,j} \left|u\_{i,j}\right|, \max\_{i,j} \left|v\_{i,j}\right|\right), 1\right) \tag{3.48}$$

The value of can be interpreted as the maximum fraction of a cell whose information can travel in one time step, multiplied by 1*.*2, which is an empirical factor that brings the approximation scheme more towards upwinding and is, from experience, more accurate [40].

With the introduction of averaged quantities

$$
\bar{u}\_{i+\frac{1}{2},j}^h = \frac{u\_{i+1,j} + u\_{i,j}}{2} \text{ and } \tag{3.49}
$$

$$
\bar{u}\_{i,j+\frac{1}{2}}^v = \frac{u\_{i,j+1} + u\_{i,j}}{2} \,, \tag{3.50}
$$

(equally for the component *v*) using the superscript *h* and *v* to indicate a horizontal or vertical averaged quantity and denoting the di↵erenced quantities likewise with *u*˜*<sup>h</sup>, u*˜*<sup>v</sup>*, etc., the nonlinear terms in equations (3.36) and (3.37) can be written as the linear combinations

$$\left(\frac{\partial u^2}{\partial x} + \frac{\partial uv}{\partial y}\right) = \frac{\partial}{\partial x}\left(\left(\bar{u}^h\right)^2 - \gamma \left|\bar{u}^h\right| \bar{u}^h\right) + \frac{\partial}{\partial y}\left(\left(\bar{u}^v \bar{v}^h\right) - \gamma \left|\bar{v}^h\right| \bar{u}^v\right) \text{ and} \tag{3.51}$$

$$\left(\frac{\partial uv}{\partial x} + \frac{\partial u^2}{\partial y}\right) = \frac{\partial}{\partial x}\left(\left(\bar{u}^v \bar{v}^h\right) - \gamma \left|\bar{u}^v\right| \bar{v}^h\right) + \frac{\partial}{\partial y}\left(\left(\bar{v}^v\right)^2 - \gamma \left|\bar{v}^v\right| \bar{v}^v\right) \,, \tag{3.52}$$

which are again approximated by centered finite di↵erence stencils.

#### 3.3.4 Boundary conditions

The boundary conditions prescribed for the model were already stated in Section 3.1.4, but since a staggered grid was introduced, applying the correct boundary conditions requires care, as some grid points lie on a boundary, while others have a boundary between them. For the no-slip boundary conditions in place, we can obtain the condition

$$u\_{i,j} = 0,\text{ if }i = 1\text{ and }n\_x + 1\tag{3.53}$$

at the west and east domain borders and

$$v\_{i,j} = 0, \text{if } j = 1 \text{ and } n\_y + 1 \tag{3.54}$$

at the south and north borders. To force zero velocity at the boundary, the conditions

$$\begin{aligned} u\_{i,j=1} &= -\
u\_{i,j=2} \quad \text{and} \\ u\_{i,j=n\_y+1} &= -\
u\_{i,j=n\_y+2} \end{aligned} \tag{3.55}$$

at the south and north domain border and

$$\begin{aligned} v\_{i=1,j} &= -\
v\_{i=2,j} \quad \text{and} \\ v\_{i=n\_x+1,j} &= -\
v\_{i=n\_x+2,j} \end{aligned} \tag{3.56}$$

at the west and east border are also prescribed. The homogenous Neumann boundary conditions for the pressure *p* yield the following conditions at the four boundaries:

$$\begin{aligned} p\_{i=1,j} &= -\ p\_{i=2,j} \\\\ p\_{i=n\_x+1,j} &= -\ p\_{i=n\_x+2,j} \\\\ p\_{i,j=1} &= -\ p\_{i,j=2} \\\\ p\_{i,j=n\_y+1} &= -\ p\_{i,j=n\_y+2} \end{aligned} \tag{3.57}$$

#### 3.3.5 Adaptations to the existing MATLAB code

Although the basic structure, which was explained in Section 3.3 and documented in [36], was used as foundation for the implemented code, much of it had to be modified in order to adapt it the problem described here.

#### External force term

As the existing code did not include outer force terms, the Boussinesq force term was added to the code and treated implicitly in equations (3.40) and (3.41).

#### Solid phase modelling

To model the liquid and solid phases in the domain, another part of the force term in the implicit equations mentioned above, was added. This can be seen as an external force density, which inhibits fluid motion in the solid part [41]:

$$\mathbf{f} = -\frac{(1-a)^2}{a^3 + 0.001} 10^8 \cdot \mathbf{u}^{\*\*} \tag{3.58}$$

37

The large number 10<sup>8</sup> was chosen to force the velocity to zero via the damping term. The liquid fraction parameter is *a* = 0 in the solid phase and rises over the mushy region to *a* = 1 in the liquid phase.

#### Phase boundaries

The boundary conditions defined in 3.3.4 are in place at the domain boundaries and the phase boundaries. To ensure this, the existing velocity field is modified to suit the conditions before the explicit time step is performed.

For the pressure *p*, this step is more complex, as the boundary conditions are implemented in the Laplace-operator to solve the Poisson equation (3.44). Therefore, the Laplace-operator must be assembled anew in every time step.

# Chapter 4 MATLAB implementation

This chapter deals with the final form of the MATLAB implementation of the numerical model described in Chapter 3. The basic structure of the code is outlined in Section 4.1, while the main functions are explained in some detail in Section 4.2. The code of these important code parts is also listed in the appendix A and referred to in this chapter.

For those interested in applying the code to conduct simulations on their own, it is highly recommended to look up Section 4.3. Therein, a detailed guide on how to specify input data for simulations is given.

It should be noted that the herein documented code was created over a period of several months. Many extensions to the initial version have been made, and the code was significantly modified several times. Although the code was written in a very comprehensible manner in the early development, measures to reduce the computational e↵ort, such as preprocessing, vectorization and the usage of the sparse matrix format, made it much more complex and probably quite dicult to read for unfamiliar users. However, it should be stressed that readability was not of crucial importance for this implementation, since this code was not developed for educational purposes, but for a designated application.

For the sake of better readability, filenames, the names of MATLAB functions or scripts and variable names are denoted by a di↵erent font: name.

## 4.1 Structure of the program

The MATLAB program developed was built as one main function with numerous sub-functions. The execution procedure is represented by the flowchart in Figure 4.1. The simulation program is executed by calling the main function PCMsolver and passing the path to the input datafile, which specifies all simulation parameters as input argument. This is discussed in detail in Section 4.3.

#### Preprocessing

In preprocessing, input data is loaded to the workspace where all important parameters are saved in the structure variable param. The finite element mesh is created by calling the function mesh2d (see 4.2.1). All data arrays are initialized, most notably, the initial temperature array T old and velocity arrays U old and V old. The function shape mat preprocesses the shape function values at the Gauss integration points. This is useful here, because all elements are of the same size and therefore the values do not change. The calculated values are stored in the structure variable gauss.

The function FEM preprocessor adds further information to this structured variable. It evaluates the values of the density, specific heat capacity and heat conductivity for every Fauss point of every element in the domain. Although the specific heat capacity must be recalculated at every time step for the PCM material, this step saves work on the time iteration. Furthermore, boundary informations are preprocessed and the relations for the sparse matrix assembly, which is carried out by the function heat2Delem new, are defined.

#### Time stepping

The time stepping is done in a for-loop, calling firstly the function heat2Delem new, which calculates the finite element matrices [*C*], [*K*] and [*R*], defined by the equations (3.24), (3.25) and (3.26). The equation system (3.23) is then solved by the function solve euler back for the nodal temperature values at the next time step T new.

If the simulation is run without considering convection, the body of the for-loop ends after saving the simulation results of the time step and reinitializing the arrays of variables T old = T new, U old = U new and V old = V new for the next time step. Otherwise, the velocities are updated as will be discussed below.

Figure 4.1: Flowchart diagram of FEM model PCMsolver

*Note.* Code given in A.1

#### Convection calculations

The velocity is updated in the function solve navierstokes, which essentially solves the Navier-Stokes equations on a staggered grid by the finite di↵erence scheme explained in Section 3.3. Although these updated velocity values are, at this point, close to zero in solid domain areas, because of the force term introduced in the implicit viscosity step (3.58), postprocessing must be performed. Therein, all velocity components, enclosed in a completely solid finite element, are set to zero. The no-slip boundary values are then newly prescribed to the updated velocities as stated by equations (3.55) and (3.56). These steps are performed for the domain boundaries and also for the liquid-solid phase boundaries, which can be observed in detail in lines 127 to 146 of listing A.1.

#### Postprocessing

After the last time step is performed, the complete simulated data is saved to the results file and optional user information is displayed.

### 4.2 Description of the key functions

In this section, the most substantial sub-functions of the main program are explained. This is done by highlighting the basic principles of implementing the model in MAT-LAB, as described in Chapter 3. The main function PCMsolver was already discussed in detail in the previous section, its listing A.1 can be found in the appendix A, together with the other key functions.

#### 4.2.1 Mesh generation in mesh2d

The finite element mesh is generated in the function mesh2d (A.2) by the call

```
[ param ] = mesh2d (param ) ;
```
in PCMsolver. As already mentioned in Section 3.2.1, the two-dimensional domain is split into an equidistant mesh of *nel* = *n<sup>x</sup> n<sup>y</sup>* four-noded bilinear rectangular finite elements, which yields the total sum of *n<sup>n</sup>* = (*n<sup>x</sup>* + 1) (*n<sup>y</sup>* + 1) nodes. These elements and nodes are arranged in the domain from left to right and from bottom to top. This is best understood by looking at Figure 4.2, which shows an example of the element and node arrangement for a mesh with *n<sup>x</sup>* = 2 and *n<sup>y</sup>* = 2. The information about the connectivity between the nodes and the elements is saved to the array param.mesh.IEN.

Figure 4.2: Image of a 2x2 mesh, with numbered nodes and elements

Furthermore, this function not only generates the coordinates of element nodes, but also those of the staggered grid points. Additionally, boundary elements are flagged, which allows quick identification when calculating the heat flux at the boundaries, and the array param.mesh.elem mat is created, which contains indexing information about the element materials based on the geometry specifications.

#### 4.2.2 Preprocessing operations

The main purpose of the functions shape mat (A.3) and FEM preprocessor (A.6) was already highlighted in 4.1. These are executed once in the simulation by the call:

```
[ gauss ] = shape mat (param ) ;
```
[ gauss ] = FEM preprocessor (T old , gauss , param ) ;

In shape mat, all occurring terms of shape matrices [*N*] and [*B*] in the system matrices (3.24), (3.25) and (3.26) are evaluated for each Gauss integration point, as well as the determinant of the Jacobian matrix (3.28). This very useful preprocessing method is possible since, in this implementation, all elements are the same size and do not change during the simulation. The shape functions are calculated by the sub-functions NmatHeat2D and BmatHeat2D, which essentially contain the formulas (3.17) and (3.15).

Function FEM preprocessor calculates further information for the FEM matrix assembly. It evaluates the values of the density, specific heat capacity and heat conductivity for every Gauss point of every element in the domain and preprocesses boundary informations, as already explained in 4.1. In addition, sparse matrix assembly relations are defined, to be used later in heat2Delem new. The sparse matrix format in MATLAB provides ecient storage of data and speeds up the processing of that data. A good example of how the element matrices are assembled to the global matrix can be found in [42], for example.

#### 4.2.3 FEM Matrix calculation in heat2Delemnew

The main FEM processing is done in heat2Delem new (A.7), where the matrices [*C*], [*K*] and [*R*], defined by the equations (3.24), (3.25), (3.26) are evaluated by an implementation of these equations. The function is called in PCMsolver by:

[C,K,R, Qin e , Qout e ,H( t ) , H l(t)] = ... heat2Delem new ( t , T old ,Uwb,Vwb, p el , gauss , param ) ;

The heat conductivity and specific (apparent) heat capacities are evaluated at each Gauss point, which is only done for elements containing PCM material. For the remaining elements, the preprocessed data of FEM preprocessor is used, which also included the preprocessed shape matrix terms. The global matrices are then assembled using the relations defined in the aforementioned function. The apparent heat capacity is calculated by the function c var (A.8), which implements definition (2.20) to account for phase change in the PCM material.

Regarding the convective term in matrix (3.25), the integration is performed by the Gauss-Lobatto quadrature instead of the common Gauss-Legendre quadrature, as already explained in Section 3.2.3. The arrays containing the velocity values *u* and *v* are in line with the scheme of the staggered grid (see Section 3.3) with already implemented boundary conditions. These values are used to interpolate the values at the Gauss points.

#### 4.2.4 Solving the matrix equations

The function solve euler back (A.10) is called thus:

```
T new = s ol v e euler back (T old ,C,K,R, param ) ;
```
It contains only one definition:

Tnew = Told + param . mesh . d t ... ( (C + param . mesh . d t K) *\*(RK Told ) ) ;

This is the MATLAB implementation of the applied backwards Euler step, yielding the formula (2.38) for the updated temperature value. The sparse matrix format allows very ecient evaluation of this step.

#### 4.2.5 Solving the Navier Stokes equations

The convection modelling discussed in detail in Section 3.3 is performed in the function solve navierstokes (A.12), which is called in the main function PCMsolver thus:

```
[ U new , V new ] = s ol v e navierstokes (U old , V old , T ij , p el , param
    ) ;
```
Therein, the force terms, which inhibit fluid motion in the solid elements, are preprocessed, as are the temperatures that are needed for the Boussinesq approximation terms. The remaining function consists primarily of an implementation of the code documented in [36].

Worth mentioning is the function Laplace p (A.13). It was developed in order to compute a finite di↵erence Laplace stencil for the pressure *p* which is consistent with homogeneous Neumann boundary conditions on the domain boundaries, as well as the internal phase boundaries.

## 4.3 Usage information

The main function PCMsolver must be run by passing a string as an input argument. This string contains the path to the MATLAB file containing the simulation input data. For instance, if the file is called input datafile.m the simulation is started by calling:

```
PCMsolver ( input datafile ) ;
```
All relevant simulation input data can be specified using such an input datafile. The code Listing 4.1 gives an example of a typical declaration structure, which is recommended. This file actually represents a typical input datafile used in the simulations carried out for the parameter study discussed in Section 6.2. In the following, both essential and optional input parameters are discussed in detail, referring to the respective lines in code Listing 4.1.

Code Listing 4.1: Listing of a typical input datafile input datafile.m

```
1 %% Input Data File 2 % The material and computational parameters are specified and written to 3 % the structure array param . 4
 5 %% name for file saving
 6 param.save.name = results ; 7
 8 %% simulate with convection 9 param.convection = 1; %(0 = off , 1 = on)
10
11 %% geometry specifications
12 % length of material domains in x - direction ( left to right )
13 param.geom.x_length_vec = [0.002 0.118];
14 % length of material domains in y - direction ( bottom to top)
15 param.geom.y_length_vec = [0.001 0.023 0.003]; %
16 % structure of simulated domain ( left to right and top to bottom )
17 param.geom.mat_mtrx = flipud([3 3
18 3 1;
19 3 3]);
20 % definition of material specifier
21 % must be in accordance with definition of material properties below 22 % PCM material specifier must be 1 23 param.geom.material = {1, pcm ; % PCM material 24 2 , stl ; % steel 25 3 , alu ; % aluminium 26 4 , air };% air 27
```

```
28 param.geom.x_length = sum(param.geom.x_length_vec); % length of PCM in x - direction 29 param.geom.y_length = sum(param.geom.y_length_vec); % length of PCM in y - direction 30
 31 %% mesh specifications
 32 nx_min = param.geom.x_length/0.0005; % minimum number of elements in x direction 33 ny_min = param.geom.y_length/0.0005; % minimum number of elements in y direction 34 nx = calc_grid(nx_min ,param.geom.x_length_vec);
 35 ny = calc_grid(ny_min ,param.geom.y_length_vec);
 36 param.mesh.nx = nx;
 37 param.mesh.ny = ny;
 38 param.mesh.nel = nx*ny; % number of elements 39 param.mesh.non = (nx+1)*(ny+1); % number of nodes 40 param.mesh.ngp = 2; % number of gauss integration points
 41 param.geom.dx = param.geom.x_length/nx;
 42 param.geom.dy = param.geom.y_length/ny;
 43
 44 %% time step specifications
 45 param.mesh.tt = 3*60*60; % total time 46 param.mesh.d_t = 0.05; % time step
 47 param.mesh.not = ceil(param.mesh.tt/param.mesh.d_t); % number of time steps
 48
 49 %% initial and boundary condition 50 % T_init [ C ] = initial temperature for element nodes 51 param.icbc.T_init =218* ones(param.mesh.non ,1);
 52 % T_in [ C ] = input temperature dependent of time 53 param.icbc.Tin = @(t,y) 230; % (t <= param . mesh .not) *200+(t > param . mesh .not) *100;
 54 % alpha_in [W / m * K] = heat conductivity dependent of y and time 55 %yh = param . geom . y_length /2;
 56 param.icbc.alpha_in = @(t,y) 700; % (y <yh) *100+(y >= yh) *20;
 57 % T_out [ C ] = output temperature dependent of time 58 param.icbc.Tout = @(t,y) 20; % output temperature [ C ]
 59 % alpha_out [ W / m *K ] = heat conductivity to environment 60 param.icbc.alpha_out = @(t,y) 1e-2; % heat conductivity [ W / m * K ]
 61
 62 %% material properties
 63 % material properties of pcm
 64 param.pcm.rho = 2050; % density [kg / m^3] ( solid )
 65 param.pcm.Tm = 220; % melting temperature [ C ]
 66 param.pcm.dTm = 1.0; % mushy region size [ C ]
 67 param.pcm.c_s = 1350; % heat capacity of solid [ J / kg* K ]
 68 param.pcm.c_l = 1492; % heat capacity of liquid [ J / kg* K ]
 69 param.pcm.lh = 108*1000; % latent heat of PCM [ J / kg]
 70 param.pcm.k_s = 0.457; % heat conductivity , solid [ W / m * K ]
 71 param.pcm.k_l = 0.435; % heat conductivity , liquid [ W / m * K ]
 72 % convective properties of pcm
 73 param.pcm.beta = 3.5e-4; % thermal expansion coeff . [1 / K ]
 74 param.pcm.grav = 9.81; % standard gravity [ m / s ^2]
 75 param.pcm.mue = 5.8e-3; % viscosity [ m^2 / s ]
 76 param.pcm.phi = 0; % angle of gravitation vector to y - axis 77 % examples of angle phi: 78 % phi = 0 ( default option ) -> gravitation in -y direction 79 % phi = pi -> gravitation in + y direction 80 % phi = pi/2 -> gravitation in -x direction 81 % phi = -pi/2 -> gravitation in + x direction 82
 83 % material properties of aluminium 84 param.alu.rho = 2700; % density [kg / m^3] ( solid )
 85 param.alu.Tm = 1000; % melting temperature [ C ]
 86 param.alu.dTm = 0; % area of melting [ C ]
 87 param.alu.c_s = 910; % heat capacity of solid [ J / kg*K ]
 88 param.alu.lh = 0; % latent heat of PCM [ J / kg]
 89 param.alu.k_s = 237; % heat conductivity , solid [ W / m * K ]
 90
 91 % material properties of steel 92 param.stl.rho = 7800; % density [kg / m^3] ( solid )
 93 param.stl.Tm = 1000; % melting temperature [ C ]
 94 param.stl.dTm = 0; % area of melting [ C ]
 95 param.stl.c_s = 540; % heat capacity of solid [ J / kg*K ]
 96 param.stl.lh = 0; % latent heat of PCM [ J / kg]
 97 param.stl.k_s = 51; % heat conductivity , solid [ W / m * K ]
 98
99 % material properties of air 100 param.air.rho = 1; % density [kg / m^3] ( solid )
```
 param.air.Tm = 1000; % melting temperature [ C ] param.air.dTm = 0; % area of melting [ C ] param.air.c\_s = 1000; % heat capacity of solid [ J / kg\* K ] param.air.lh = 0; % latent heat of PCM [ J / kg] param.air.k\_s = 0.024; % heat conductivity , solid [ W / m \* K ] 106 %% data recovery param.save.intervall = 6\*ceil(1/param.mesh.d\_t); % time step intervall <sup>109</sup> % for T ,U , V to be saved <sup>110</sup> param.save.make\_snapshot = 1; %(0 <sup>=</sup> off , 1 = on) % dump snapshots of data <sup>111</sup> param.save.snapshot\_int = 120; % time step intervall to make snapshots 112 %% display information for user <sup>114</sup> % (0 <sup>=</sup> off , 1 = on) param.display.waitbar = 1; % show waitbar during run <sup>116</sup> % info : each waitbar update needs about 0.2 seconds computation time <sup>117</sup> param.display.update = 10; % minimum time between waitbar updates in seconds <sup>118</sup> param.display.Tplot = 0; % plot temperature distribution at the end <sup>119</sup> param.display.Vplot = 0; % plot velocity distribution at the end <sup>120</sup> param.display.matplot = 0; % plot visualisation of material in domain

As mentioned above, all input parameters are specified in the structure array param. For more clarity, the fields of param again contain structure arrays.

In line 6, the basic string for the filenames of the result files is declared. Line 9 contains binary information about whether the e↵ect of natural convection is to be considered in this simulation.

#### Geometry specifications

The geometry specifications are defined in lines 11 to 30. The length of di↵erent material domains is specified using vectors and the organisation can be easily understood by looking at Figure 6.1. The field param.geom.mat mtrx contains a matrix of indices specifying certain materials. These indices are to be defined arbitrarily in param.geom.material, except for the fixed index 1, which refers to the PCM material. The syntax of the matrix formulated material declaration chosen here can again be best put into context by observing Figure 6.1.

#### Mesh specifications

The essential mesh specifications are defined in param.mesh.nx and param.mesh.ny (lines 36 and 37), which are the number of elements in x direction and y direction respectively. Lines 32 to 35 may be used by calling the created function calc grid to calculate a number of elements that is compatible with the material domains specified in lines 13 and 15, in the means that no element contains more than one material.

#### Time step specifications

Lines 45 and 46 contain the time step specifications, which are the total simulation time param.mesh.tt and the time step size param.mesh.d t.

#### Initial and boundary conditions

The possible values for initial and boundary conditions are defined in lines 49 to 60, which are the heat conductivities ↵*in* and ↵*out* and temperatures *Tin* and *Tout* belonging to the Robin boundary conditions on the west and east side of the domain. The initial temperature field is prescribed for each element through the vector param.icbc.T init.

#### Material properties

Material properties are defined in structure arrays and saved in fields of param, of which the names must be declared in accordance with the definition of material specifiers in param.geom.mat mtrx. The basic thermophysical material properties for the phase change material are defined in lines 64 to 71, and so is also the mushy region size ✏ in line 66. The same can be done for other materials, as given here for aluminium, steel and air in lines 84 to 105. The convective simulation parameters are defined in lines 73 to 76 and contain the thermal expansion coecient param.pcm.beta, the viscosity param.pcm.mue and the angle param.pcm.phi. This angle, which is defined as the angle of the gravitational vector to the y-axis, was introduced for the parameter study in Section 6.2 to study di↵erent domain orientations. For the intuitive case that the y-axis is aligned to the gravitational force vector, this angle takes the value zero. Other values can be defined in radian units in a clockwise manner, as given in the example values in lines 78 to 81.

#### Data recovery

Some options are available regarding the recovery of the simulation data of the node temperature values and velocity value matrices. param.save.intervall defines the number of time steps between moments at which the mentioned quantities are to be saved. This is useful to keep the amount of aggregated data moderate. Snapshots of the data can be saved by setting the binary variable param.save.make snapshot to the value 1, at intervals defined by the associated variable param.save.snapshot int. This option allows insights into the results of ongoing simulations.

#### Display parameters

With the binary variable param.display.waitbar, the display of a waitbar, showing the progress of the simulation and remaining time until completion, can be switched on or o↵. The corresponding variable param.display.update declares the time between waitbar updates. However, this option comes with the disadvantage of significant computational e↵ort for small update intervals, because of the rather slow MATLABinternal function draw. More visualization options are available in lines 118 to 120, which can be enabled by setting them to the value 1 and disabled by setting them to the value 0.

# Chapter 5

## Validation

In this chapter, the results of the model validation for the latent heat storage model are presented. The validation was conducted in three steps.

The first validation step was to compare the results of a simulation where only heat conduction was considered to those of the analytical solution of the Stefan problem, introduced in Section 2.1.4.

The second test case is the classic problem described in [5], that is phase change in a cavity filled with gallium with buoyancy forces. This test case aimed at validating the convection model described in Section 3.3.

In the next step, a cross-validation with the CFD software ANSYS Fluent was performed. As part of the parameter study, which will be discussed in Chapter 6, the model results are compared to those replicated with the Fluent software package.

Another test case was performed which was lass a validation attempt than a test of the model stability and feasibility in certain simulation scenarios. This test scenario alternated between charging and discharging a PCM storage to see if numerical problems arose and if the results complied with physical principles.

For the two validation attempts of the heat conduction model and convection model, a mesh convergence analysis was performed to study the e↵ect of the discretization parameters mesh size and time step on the numerical solution. Also, the e↵ect of the chosen mushy region temperature range and possible complications are scrutinized, as are the e↵ect of the aspect ratio of the chosen finite elements.

### 5.1 Validation of heat conduction model

To validate the accuracy and convergence of the model without convection, the onedimensional Stefan problem was considered, as described in Section 2.1.4. The material parameters used in this simulation are those of KNO3-NaNO<sup>3</sup> given in Table 3.1, but to simplify the analytical solution, the specific heat capacity *c<sup>l</sup>* = *c<sup>s</sup>* and thermal conductivity *k<sup>l</sup>* = *k<sup>s</sup>* were set to the values of the solid substance. The melting point was defined as *T<sup>m</sup>* = 220 *C* and the left wall input temperature as 235 *C*. The initial temperature of the simulated domain was naturally set to *Tinit* = *T<sup>m</sup>* ✏, with the half width of the mushy region ✏. The heat conduction was simulated in a domain with the length of 50 *mm* and over a period of five hours, 5 *h*, in which the melting front could not quite spread over the whole domain.

This validation simulation was carried out for di↵erent values of the time step size *t*, mesh size *x* and mushy region parameter ✏, where in each case the other two parameters were held at a constant value. Therefore, the influence of the mentioned numbers on the result was studied separately as given in the sections below.

The RMSE (root mean square error) statistical values [43]

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left( s\_{analytic} - s\_{simulation} \right)^2} \tag{5.1}$$

were used to compare analytical predictions and simulation results, where *s* is the position of the melting front and *n* the number of comparison points. The melting front positions were evaluated and compared at *n* = 1800 points in time, spread equidistantly over the five hours of simulation time.

To quantify the convergence of the solution with finer mesh and smaller time steps, the melting front positions were also compared to the most precise simulation in that regard. For this purpose, the values *sanalytical* in equation (5.1) are compared with those of the most precise simulation.

#### E↵ect of mesh refinement

To study the e↵ect of mesh refinement, the Stefan problem was simulated for di↵erent values of the mesh size *x*, while the time step size *t* = 1 *s* and mushy region parameter ✏ = 1 *C* were held constant. The RMSE compared to the analytical prediction and the RMSE compared to the simulation with *x* = 0*.*1 *mm* are given in Table 5.1.

Table 5.1: RMSE of melting front position for di↵erent values of mesh size or number of elements compared to Stefan solution (analytical) and *x* = 0*.*1 *mm* simulation (convergence)


With the RMSE compared to the most precise simulation with *x* = 0*.*1 *mm*, we can observe very fast convergence with decreasing mesh size. However, it can be noted that, compared to the analytical solution, the simulation does not become more accurate with decreasing mesh size after *x* = 2 *mm*. Obviously, the accuracy limiting factor, in this case, lies in other simulation parameters.

#### E↵ect of applied time step

In the next step, the mesh size *x* = 0*.*5 *mm* and mushy region parameter ✏ = 1 *C* were held constant, while the simulation of the Stefan problem was carried out for di↵erent values of the time step size *t*. Again, the RMSE values are given in Table 5.2, in the same fashion as in Section 5.1.


Table 5.2: RMSE of melting front position for di↵erent values of time step size compared to Stefan solution (analytical) and *t* = 0*.*1 *s* simulation (convergence)

Again, the RMSE compared to the most precise simulation with *t* = 0*.*1 *s* indicates very fast convergence with decreasing time step size. The error compared to the analytical solution remains at a rather steady value for time steps *t* = 5 *s* and smaller. This, however, can be very di↵erent for other mushy region sizes ✏. For smaller values of the latter, the time step size should also be decreased. Otherwise, the phase change is partly or completely skipped in certain finite elements, which is a known problem of the apparent heat capacity method and which will be clearer in the next subsection.

#### E↵ect of mushy region size

The e↵ect of the mushy region size ✏ on the simulation result of the Stefan problem was examined under constant values of mesh size *x* = 0*.*5 *mm* and time step size *t* = 1 *s*. The error values of the melting front position compared to the analytical solution are given in Table 5.3.

Table 5.3: RMSE of melting front position for di↵erent values of mushy region size ✏ compared to Stefan solution (analytical)


The quantities in Table 5.3 are best put into context by looking at the time evolution of the melting fronts of these di↵erent simulations, given in Figure 5.1. The melting front lags slightly behind the analytical solution for large mushy region size values, which can be explained by the distinctly altered actual physical problem, the Stefan problem, where no mushy region exists. Having said this, for small values of the mushy region size, the simulation does not become accurate, due to the nature of the apparent heat capacity method. As mentioned above, the phase change is partly or completely skipped in certain finite elements if the mushy region is small and the time step is large. This causes a higher RMSE value for the simulations with mushy region size ✏ = 0*.*2 *C* and ✏ = 0*.*1 *C*, compared to the simulation with ✏ = 0*.*5 *C*. In addition, the e↵ect of this error can be seen very clearly for the simulated melting front with mushy region size ✏ = 0*.*1 *C*, which is plotted as a dashed red line in Figure 5.1.

Figure 5.1: Time evolution of the melting front of simulations with di↵erent values of the mushy region size ✏ compared to the melting front of the analytical solution

Looking at Figure 5.1, we can obtain that the error values for the simulations with mushy region size ✏ = 1 *C* and ✏ = 0*.*5 *C* given in Table 5.3 should be taken into consideration with care. These small deviations from the analytical solution most probably arise from two competing sources of error, which cancel each other out. Firstly, phase change is partly skipped in some elements during the early stage of the simulation with high temperature gradient, which causes the melting front to be slightly more advanced, compared to the analytical solution. Secondly, the physical problem is distinctly altered by introducing a mushy region, which causes a slower progress of the melting front, which therefore cancels out the first mentioned error.

Based on the findings of the last section, it is suggested that the magnitudes of the time step and mushy region size be chosen so that the mushy region always overlaps its previous position, ensuring that no finite element skips the phase change. Due to the nature of the apparent heat capacity method, these elements would not be exposed to the e↵ect of latent heat. These parameters should thus be chosen with care, and the findings of this section might help to give a basic idea of the necessary scaling.

#### Final 1D validation result

A final 1D simulation was carried out, aiming for an excellent prediction of the Stefan problem. For this, the mesh size was chosen as *x* = 0*.*5 *mm*, which was already seen as a good approximation. The time step size *t* = 0*.*02 *s* was chosen very small to easily handle a small mushy region of ✏ = 0*.*1 *C*.

The RMSE value compared to the analytical solution was found to be 1*.*08 10<sup>4</sup> *mm* but the absolute deviation decreased with growing simulation time. Thus, the relative error of the melting front prediction was well under 1 %, except for the early stages of the simulation. On the basis of the above, it can be concluded that the developed model gives an excellent approximation of heat conduction problems if the simulation parameters are carefully chosen.

## 5.2 Validation of convection model

#### 5.2.1 Test case

The validation of the convection model is based on the report [5], in which the role of natural convection in solid-liquid-interface motion is studied during melting and solidification of the metal gallium on a vertical wall.

The rectangular test cell (8.89 cm in height, 6.35 cm in width) was heated on one side and a heat sink was located on the opposite side wall. The other walls were insulated with Plexiglass, and the front and back sidewalls contained an air gap between two Plexiglass plates. These sidewalls were assumed to be adiabatic in the validation simulations. In the experiment used for this validation, the heat sink was held at a constant temperature of 28*.*3 *C*, which is also the initial temperature of the gallium in the domain, whereas the heated wall was held at constant 38*.*0 *C*. Since the mathematical model of this problem is naturally described with Dirichlet boundary conditions, prescribing a constant temperature at the walls, but Robin boundary conditions are implemented in the FEM model, the thermal conductivity at the wall was chosen to have a very high value ↵*boundary* = 10<sup>10</sup> W/mK.

The gallium used in the experiments [5] had a purity of 99*.*6 % and a melting temperature of 29*.*78 *C*. All material properties needed for the simulation of this test case were extracted from the article [44], where a numerical investigation of exactly this experiment was conducted. A total of 26 equally spaced thermocouples were installed on the top and bottom walls and another 17 thermocouples along the center line, to measure the horizontal temperature distributions and therefore allow for tracing the interface of the melting front.

The shape of the melting front is given at nine di↵erent moments in the experiment, namely at minutes 2, 3, 6, 8, 10, 12.5, 15, 17 and 19 and was extracted from the graphic in the paper. From these melting fronts, the melted fraction in the gallium cavity was deduced and used as a quantity of comparison. Since the shape of the melting front, calculated from simulation results, is not a mesh-independent solution, the liquid fraction is calculated via the current latent heat stored in the cavity. To draw comparisons, again the RMSE (root mean square error) value

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left( a\_{experiment} - a\_{simulations} \right)^{2}} \tag{5.2}$$

is used with the *n* = 9 time comparison points, where *a* is the liquid fraction in the cavity.

#### 5.2.2 Validation and convergence analysis

A convergence analysis was performed on the e↵ect of mesh resolution, timestep size, mushy region temperature range and also aspect ratio on the simulation accuracy compared to the experimental solution in the explained test case (see Section 5.2.1). Thus, a validation for di↵erent values of these parameters is given.

#### E↵ect of mesh refinement

To study the e↵ect of mesh refinement, the gallium test case, as described above, was simulated for di↵erent values of the mesh size, while the time step size *t* = 25 *ms* and mushy region parameter ✏ = 0*.*1 *C* were held constant. The mesh was chosen so that *x*=*y*. Therefore, the particular mesh size is labelled only by the value of *x* in the following. The RMSE compared to the experimental data and the RMSE compared to the simulation with *x* = 0*.*5 *mm* are given in Table 5.4.

With the RMSE compared to the most precise simulation with *x* = 0*.*5 *mm*, we can observe a significant improvement from the coarsest mesh to the mesh with *x* = 2*.*0 *mm*. From this value on, the convergence is very slow and the simulations do not considerably increase in accuracy compared to the experimental data.

Table 5.4: RMSE of liquid fraction for di↵erent values of mesh size or number of elements compared to test case (experimental) and *x* = 0*.*5 *mm* simulation (convergence)


#### E↵ect of applied time step

The gallium test case was simulated for di↵erent time step sizes, while the mesh size *x* = *y* = 1 *mm* and mushy region temperature range ✏ = 0*.*1 *C* were held constant. The RMSE compared to the experimental data and the RMSE compared to the finest simulation with *t* = 2 *ms* are shown in Table 5.5.

Table 5.5: RMSE of liquid fraction for di↵erent values of the time step size compared to test case (experimental) and *t* = 2 *ms* simulation (convergence)


Looking at the RMSE values, convergence of the simulation result with decreasing time step size was found, although not quite as fast as with refinement of the mesh. The same can be said about the improvement of accuracy compared to the experimental data. Still, a root mean square error of below 2 %, which accounts to a relative deviation of around 3 % of the experimentally obtained liquid fraction, can be seen as a very good validation of the convection model.

#### E↵ect of mushy region size

The e↵ect of the mushy region size ✏ on the simulation result of the gallium test case was examined under constant values of mesh size *x* = 1 *mm* and time step size *t* = 25 *ms*. Again, the RMSE compared to the experimental data is shown in Table 5.6. The RMSE values in Table 5.6 are similar to the findings in Section 5.1. Again, the melting front lags slightly behind the compared solution, which in this case is the liquid fraction, for large values of the mushy region size, but it is not causing large errors. For small values of the mushy region size, the simulation gets poorer, for the same reasons already explained for the 1D validation. We can therefore conclude that when using the apparent heat capacity method for isothermal phase change problems, approximated by a very small mushy region, a very small time step size must also be used.

Table 5.6: RMSE of liquid fraction for di↵erent values of the mushy region size compared to test case (experimental)


#### E↵ect of the aspect ratio

Another question with respect to stability and convergence of the developed model was whether large aspect ratios *<sup>x</sup> <sup>y</sup>* or *<sup>y</sup> <sup>x</sup>* cause any problems for the simulation. Such a choice of mesh might be convenient for simulation domains with large aspect ratios, because of the reduced total number of elements. To answer this question, the gallium test case was simulated for di↵erent values of the aspect ratio, while the time step size *t* = 25 *ms* and mushy region parameter ✏ = 0*.*1 *C* were held constant. The melted volume fraction of such simulations, where the aspect ratio was tested with respect to larger values of *x*, is shown in Figure 5.2. The same analysis was done for larger values of *y*, and the melted volume fraction is shown in figure 5.3.

Figure 5.2: Time evolution of liquid fraction for di↵erent aspect ratio values *ar* <sup>=</sup> *<sup>x</sup> y* relative to the base case *x* = *y* = 0*.*5 *mm* compared to experimental data with variable *x*

It can be seen that the deviation from the experiment gets larger with increasing aspect ratio, but still does not exceed the range of error of the coarsest mesh *x* = 4 *mm*, *y* = 4 *mm*, which is denoted by the aspect ratio 8/8 in both figures mentioned. Therefore, this approximation error can be traced back to the coarseness of the mesh, and no complications that come with large aspect ratios have been found.

Comparing Figures 5.2 and 5.3, we can obtain that the refinement of *x* has the same e↵ect as the refinement of *y*. Furthermore, it can be observed that reasonably fine meshes already produce good simulation results, which was previously stated in Section 5.2.2.

It should be noted though, that it cannot be ruled out that complications during simulations with large aspect ratios of the mesh could occur for di↵erent domain designs, and it is therefore recommended to test this in an analysis of mesh convergence, as was conducted in this chapter.

Figure 5.3: Time evolution of liquid fraction for di↵erent aspect ratio values *ar* <sup>=</sup> *<sup>x</sup> y* relative to the base case *x* = *y* = 0*.*5 *mm* compared to experimental data with variable *y*

#### Final validation results of the convection model

To conclude the validation of the convection model via the gallium test case, the shape of the melting front for selected simulations is discussed here.

In Figure 5.5, the melting fronts of two of the simulations already described in Section 5.2.2, with di↵erent mesh sizes *x*=*y*, a time step size of *t* = 0*.*25 *ms* and mushy region parameter ✏ = 0*.*1 *C* are plotted at di↵erent comparison moments. Although the RMSE compared to the experimental data of these two simulations (see Table 5.4) was similar, there is still a significant di↵erence in the shape of the developing melting fronts. The same is true for the two simulations with *t* = 5 *ms* and *t* = 50 *ms*, both with *x*=*y* = 1 *mm* and ✏ = 0*.*1 *C*, which are given in Figure 5.4. Interestingly, the shape of the melting front of the coarser simulation seems to fit the experimentally obtained melting front better than the one with the finer parameters. A possible explanation of this is the separation of physical and numerical e↵ects, whereby expected modelling and discretization errors might cancel each other out on a coarse mesh compared to a finer mesh solution [45]. This e↵ect of course is already well known in CFD simulations and was reported numerous times, in the reference [46], for example, where it was found that reasonably coarse meshes might yield better results compared to finer mesh structures and that further refinements might even deteriorate the results.

Figure 5.4: Melting fronts of two simulations with di↵erent time step size compared to experimental data, given at five time points

In spite of the peculiarities explained in this section, overall a good agreement of the model predictions with the experimental results published in [5] was found.

Figure 5.5: Melting fronts of two simulations with di↵erent mesh size compared to experimental data, given at five time points

### 5.3 Cross-validation with ANSYS Fluent

In the course of conducting the parameter study, which will be discussed in Chapter 6, simulations were carried out in the developed MATLAB model and also run in ANSYS Fluent. Thus, a cross-validation with the renowned CFD software is performed, a brief illustration of which is shown here.

In the cases which could be compared using the two simulation methods, a maximum relative deviation of 13 % at any time step in the simulated liquid fraction was found when compared to the prediction of Fluent. In general though, the observed discrepancy was considerably smaller than that. As an example, the comparison for the case of a standard PCM cavity (as will be described in Chapter 6) in horizontal position and heated from the left side, shall be given here. Figure 5.6 shows the phase structure in the cavity at di↵erent times calculated by Fluent on the left side and calculated by the MATLAB model on the right side. Not only the extent of the liquid phase, but also the shape of the melting front, simulated by the developed FEM model, were found to be in very good agreement with the results from Fluent.


(g) after 180 minutes

Figure 5.6: Phase structure in PCM cavity at di↵erent times

*Note.* The solid region is portrayed in blue and the liquid region is portrayed in red for the Fluent simulation (left side) and yellow for the FEM simulation (right side).

## 5.4 Test of a fluctuating charging/discharging scenario

To determine the stability of the developed model in a simulation of operating latent heat hybrid storage, a test scenario with multiple irregular charging and discharging cycles was performed. The goal was to discern whether two or even more separated liquid cells could develop in this PCM storage and to observe whether any numerical problems thereby occur.

It should be pointed out that, the CFL condition, which is introduced when the nonlinear Navier-Stokes term is treated explicitly (see Section 3.3.2), is necessary for the stability of the simulation. This condition implies that the maximum velocity occurring in the domain must be smaller than the ratio of the spatial resolution to the time step size.

The geometry and material properties used in this examination are those of the gallium validation case, described in Section 5.2. The input temperature on the left side wall of the domain was specified di↵erently, and in an irregularly fluctuating manner, which is shown as a dashed line in Figure 5.7. The mean temperature in the gallium cavity, as well as the maximum velocity appearing in the molten gallium, are presented in this Figure as solid blue and orange lines respectively.

It can be found that the velocity in the cavity builds up very quickly to the maximum value and is later diminished by declining driving forces, namely, temperature di↵erences in the domain. This reaction can be seen for examples at around 7, 27 and 37 minutes, when the gallium is fully liquefied, and when the input temperature drops and a solid layer develops at the left side wall. Hence, no further adaption to the model to slow down the convective motion must be made. Furthermore, no numerical diculties or 'unphysical behaviour' was observed during this test simulation scenario.

As to the development of separated liquid cells in the gallium cavity, no more than two of them could be observed at any time in the simulation. The dynamic of the convective heat transfer leads to the separated zones merging very quickly and therefore makes more than two of these zones coexisting rather unlikely. Although these improbable scenarios might occur in very wide and flat cavities, the shape of the melting front studied later in Chapter 6 leads to the assumption that they will not occur in this kind of cavity with aluminium fins. An example snapshot of the temperature and velocity field during the merge of two separated liquid cells is shown in Figures 5.8 and 5.9. The solid fraction can be easily obtained from the representation of the velocity distribution in 5.9.

Figure 5.7: Temperature and velocity evolution over time in fluctuating gallium test case charging/discharging scenario

Figure 5.8: Temperature distribution in the gallium cavity at the simulation time of 26 minutes during charging/discharging scenario

Figure 5.9: Velocity distribution in the gallium cavity at the simulation time of 26 minutes during charging/discharging scenario

## Chapter 6

## Parameter studies

As an example of what can be done with the developed FEM MATLAB model, a parameter study conducted to investigate di↵erent designs specifics of latent heat storage cavities is shown here. Also, the dynamics of melting and solidification for di↵erent orientations of a selected PCM cavity were examined for the charging and discharging mode. For this selected cavity design, the e↵ect of natural convection is presented in contrast to simulations which do not consider convection.

## 6.1 Study of di↵erent cavity designs

As already mentioned in this thesis, a number of di↵erent designs are currently under consideration for the PCM cavity of the HyStEPs hybrid latent heat storage. At this point, the properties of the materials and operation parameters are not yet known. The parameter study of di↵erent cavity dimensions discussed in this section was conducted at a very early stage in the project. Many properties were still unknown or open for discussion. Also, natural convection was not yet implemented into the model, which is why this parameter study considers only heat conduction. Therefore, the simulations conducted and the results obtained might not be of great value for future considerations. However, some unique characteristics of the cavity design discussed here could also be found for di↵erent dimensions of the cavity simulated here. Furthermore, the methodical approach to accomplish first estimations of the behaviour of chosen designs can be adopted in the future.

The basic design of the latent heat storage in consideration is that of an aluminiumframed PCM cavity attached to the steel containment of the steam storage, while the outward-facing side of the cavity is insulated. Of course, the real life implementation will most probably look di↵erent due to engineering reasons, but nevertheless, should be well approximated by this model. To enhance the heat transfer in the PCM cavity, aluminium fins should be placed in the plane of the axis of the cylindric steam storage. A sketch of such a fin section is shown in Figure 6.1.

Vienna University of Technology Dominik Pernsteiner 1 Figure 6.1: Sketch of a typical PCM cavity fin section attached to the steel containment of the steam storage

#### 6.1.1 Simulated parameters

The simulation domain of the PCM cavity was defined as shown in Figure 6.1, with sections of di↵erent materials. The material properties are those shown in Table 3.1. The width of the steel containment was assumed to be 10 *mm*, and the heat transfer from the steam to the steel containment was described by the thermal conductivity ↵*in* = 100 *<sup>W</sup> <sup>m</sup>*2*<sup>K</sup>* . The heat transfer through the well insulated outwards facing wall was described by the thermal conductivity ↵*out* = 0*.*02 *<sup>W</sup> <sup>m</sup>*2*<sup>K</sup>* . Again, these parameters might be very di↵erent in reality and should be thoroughly discussed for future evaluations.

Regarding the computational parameters for these studies, a mesh of 80 times 80 elements was chosen. The mushy region parameter was set to ✏ = 2*.*0 *C*, together with a time step size of *t* = 10 *s*. The initial temperature in the domain was set to 190 *C* and the input temperature was set to 235 *C* for one hour of time and back to 190 *C* for a second hour, which should represent a simplification of a typical operation cycle of the steam storage, where melting and solidification can be observed.

The aim of the conducted parameter study was to obtain an 'ideal' value of the geometric aspects of a fin section, meaning the dimensions of the PCM *XPCM* and *YPCM* and the width of the aluminium fins *Yaluminium* for this kind of storage operation. The width of the aluminium between the steel containment and the PCM was set to a constant *Xaluminium* = 5 *mm* and the steel containment width was set to a constant *Xsteel* = 10 *mm*.

A total of 32 simulations were carried out, of which 16 simulated a fin width of *Yaluminium* = 5 *mm* and 16 a fin width of *Yaluminium* = 2*.*5 *mm*. The dimensions of the PCM were varied between *XPCM, YPCM* = 5 *mm* and *XPCM, YPCM* = 50 *mm*. The simulated cases together with the associated tagging are given in Tables 6.1 and 6.2.


Table 6.1: Simulated parameter study cases and associated tags with fin width *Yaluminium* = 5 *mm*


Table 6.2: Simulated parameter study cases and associated tags with fin width *Yaluminium* = 2*.*5 *mm*

#### 6.1.2 Simulation results

To analyse the simulation results and compare the di↵erent cavity dimensions, the evolution of the enthalpy change per surface area *<sup>H</sup> <sup>A</sup>* was calculated, where *<sup>A</sup>* is the surface area between the PCM fin section and the steam storage. The surface therefore takes the value

$$A = \left(Y\_{PCM} + 2Y\_{aluminum}\right) \cdot 1m\text{ .}\tag{6.1}$$

In Figures 6.2 and 6.3, this quantity and the latent heat fraction of it is plotted for the di↵erent simulated cases. In the parameter studies conducted here, the latent heat fraction of the total enthalpy change *H* is evidently relatively small. This is because the containment of 10 *mm* steel and 5 *mm* aluminium was also considered, which o↵ers a larger storage of sensible heat.

It can be obtained that, for the selection of longer cavities (*XPCM* = 50 *mm,* 25 *mm*), the lower cavities can store more enthalpy per surface area and the charging and discharging process is also quicker. This di↵erence in speed of charging and discharging between the cavity heights can also be obtained for the shorter cavities, as long as

(c) cases with constant *XPCM* = 10 *mm*

(d) cases with constant *XPCM* = 5 *mm*

## Figure 6.2: Enthalpy per surface area of the simulated cases with fin width *Yaluminium* = 5 *mm*

*Note.* The total enthalpy change *H* is plotted using solid lines, the latent heat fraction using dashed lines.

the PCM is not fully liquefied. At that point, the higher cavities are able to store more enthalpy because of the higher ratio of PCM to aluminium. After all, the dimensioning of the cavity is also a question of cost. If unused PCM material must be

(a) cases with constant *XPCM* = 50 *mm*

(b) cases with constant *XPCM* = 25 *mm*

(c) cases with constant *XPCM* = 10 *mm* (d) cases with constant *XPCM* = 5 *mm*

Figure 6.3: Enthalpy per surface area of the simulated cases with fin width *Yaluminium* = 2*.*5 *mm*

*Note.* The total enthalpy change *H* is plotted using solid lines, the latent heat fraction using dashed lines.

avoided because of budgetary reasons, a cavity as in case 41v10 should be preferred. If the material costs are not considered and the total storable enthalpy should be maximized, cavity dimensions as in case 11v40 are the optimal choice.

Looking at Figure 6.4, we can conclude that the smaller fin width of *Yaluminium* = 2*.*5 *mm* should be preferred over the version with *Yaluminium* = 5*.*0 *mm* because the over all storage capacity is higher and the charging/discharging speed is not diminished by the thin fins. This could of course be di↵erent for PCM materials with higher thermal conductivity, in which case the heat transfer could be limited by the width of the aluminium fins. However, convection, which leads to an e↵ectively higher heat conductivity, could also provoke this e↵ect. Finally, constructive constraints must also be considered, which generally include a lower limit for the minimum width of aluminium fins.

(a) cases with constant *XPCM* = 50 *mm* and *YPCM* = 5 *mm* (b) cases with constant *XPCM* = 5 *mm* and *YPCM* = 50 *mm*

Figure 6.4: Enthalpy per surface area of the simulated cases with di↵erent fin width *Yaluminium*

*Note.* The total enthalpy change *H* is plotted using solid lines, the latent heat fraction using dashed lines.

Considering the width of the fins, it should be stressed that the results obtained in this study should not lead to the assumption that thinner fins are always the better choice. In fact, further investigations on the fin dimensions of a cavity, as described in Section 6.2.2 with *XPCM* = 118 *mm*, showed a significant influence of the fin width on charging/discharging speed. Although these investigations are not particularized in this thesis, this should again emphasize the importance of further studies on the dimensioning of the HyStEPs cavity.

## 6.2 Comparison of di↵erent cavity orientations

In this section, the e↵ect of the orientation of a typical PCM cavity with respect to its position on the steam storage is studied. The term 'typical' PCM cavity is of course very subjective, in that it is problem-dependent and based on many di↵erent assumptions. As mentioned earlier, the cavity dimensions studied in Section 6.1 were later found to be of little relevance for future application. Based on an assessment of fellow HyStEPs project partners, the amount of PCM needed to produce a desired latent heat storage capacity for a designated application case could be estimated and leads to a required PCM coating of the steam storage of around 120 *mm* thickness. Thus, the geometry of an, at this point, more representative PCM cavity fin section, described in detail in the following, is simulated in this section.

#### 6.2.1 Simulations in ANSYS Fluent

The simulations carried out here in the developed MATLAB model were also run in ANSYS Fluent. This was done to double-check the results and to have a second tool ready and parametrized for future simulations. In Section 5.3, the successful crossvalidation was already reported. As the simulations in Fluent are not essentially part of this thesis, these will not be discussed in detail. It should only be noted, that all results obtained from the two simulation methods were very similar (see for example 5.6) and can be interpreted the same way.

Another observation worth mentioning is the negligible e↵ect of small mushy regions for the kind of phase change problems simulated. A comparison of a typical parameter study case, as reviewed in this section and simulated in Fluent once with a mushy region parameter of 0*.*1 *C* and once with isothermal phase change (which Fluent can handle), showed a relative di↵erence in the results of less than 1*.*5 %. This is a further convenience of the developed MATLAB model with regard to the apparent heat capacity method, because it means that isothermal materials can also be approximated fairly well by a very small mushy region.

#### 6.2.2 Simulation parameters

The simulated domain still represents the PCM cavity fin section illustrated in Figure 6.1, except for the steel containment, which was not implemented in this study in order to lower the computational e↵ort. This can be justified by the finding during the preliminary parameter study, explained in Section 6.1, that the steel containment has a negligible e↵ect on the transient behaviour of the whole cavity. The dimensions of the fin section took fixed values in the simulations described here. The PCM was defined by *XPCM* = 118 *mm* and *YPCM* = 23 *mm*, the width of the aluminium fins was set to *Yaluminium* = 1 *mm*, and between the PCM and the steam storage an aluminium layer of *Xaluminium* = 2 *mm* was defined.

The material properties used for these simulations can again be found in Table 3.1. The heat transfer from the steam to the aluminium layer (via the neglected steel containment) was described by the thermal conductivity ↵*in* = 700 *<sup>W</sup> <sup>m</sup>*2*<sup>K</sup>* . The heat transfer through the outward facing wall was assumed as adiabatic.

Regarding the computational parameters, a convergence analysis with respect to mesh and time step size was performed. For the chosen mushy region parameter of ✏ = 0*.*1 *C*, convergence was reached with a time step size of *t* = 0*.*1 *s*. A mesh with quadratic elements of 0*.*5 *mm* length was chosen for the simulations, which agreed best with the comparative results from ANSYS Fluent.

For the simulation of the charging process of the PCM fin section, the initial temperature in the domain was set to 218 *C* and the input temperature was set to 230 *C* for three hours of simulation time. For the discharging simulation, an input temperature of 200 *C* was assumed and the initial temperature set to 230 *C*.

HyStEPs: Location of a Fin Section

Vienna University of Technology Dominik Pernsteiner 21 Figure 6.5: Schematic illustration of the five simulated cavity orientations in relation to the gravity vector

The simulations for charging and discharging of the PCM cavity, which included both the modelling of heat conduction and natural convection, were carried out for five di↵erent positions on the steam storage with characteristic angles in relation to the gravity vector. The five di↵erent orientations 0, 45, 90, 135 and 180 are illustrated in Figure 6.5.

#### 6.2.3 Simulation results

The change of total enthalpy during the charging process, as well as the latent heat fraction, are presented in Figure 6.6. Another representation of the simulation results is shown in Figure 6.7, where the evolution of the melted volume fraction during the charging and discharging process of the cavity is charted for the di↵erent orientations, as is the case in which only heat conduction is considered.

The horizontal cavity (90) just reaches total liquefaction and, therefore, the fully charged state, only at the very end of the three hour simulation time, while for the upward facing cavities (0, 45), this state is already reached at two and a half hours. The downward facing cavities (135 and 180) could not be completely charged during the three hours of simulation time and a solid PCM fraction of more than 10 % remained. Furthermore, these two downward facing cavities show almost the exact same transient behaviour and only slightly outperform the charging process of the simulated case, where natural convection was neglected.

Figure 6.6: Change of enthalpy during the charging process

*Note.* The total enthalpy change *H* is plotted using solid lines, the latent heat fraction using dashed lines.

Looking at Figure 6.7(b), it can be concluded that the orientation had no influence on the performance of the discharging process and that the latter is e↵ectively the same as for the case in which natural convection was neglected.

In the appendix B, some plots of the temperature and velocity distribution of the simulation results for the five di↵erent cavity orientations are given, to illustrate the discussion in this chapter.

Figure 6.7: Evolution of the melted volume fraction during charging (a) and discharging (b) of the cavity for di↵erent orientations

## 6.3 Discussion of conduction vs. convection simulations

Based on the results given in the previous section, we can conclude that the e↵ect of natural convection on the charging process of the PCM cavity is highly dependent on the respective orientation. It cannot be neglected for upward facing cavities and is still a significant factor for cavities in horizontal position. For downward facing cavities however, the e↵ect of natural convection is negligible and its behaviour can be described by heat conduction only, although it may be possible to improve accuracy by introducing a convective enhancement factor [47].

During discharging, the e↵ect of natural convection was found to be insignificant, regardless of the cavity orientation.

# Chapter 7 Conclusion

This last chapter presents a short summary of the overall approach to the stated problems and the main findings during the work on this diploma thesis. The findings are evaluated and put into perspective using previous research and possible applications. Finally, future investigation objectives and extensions to the developed MATLAB model are suggested.

## 7.1 Summary

The main task of this diploma thesis was the development of a MATLAB model, preferably via the finite element method, to solve phase change problems in typical latent heat storage geometries, considering both heat conduction and natural convection. To this end, comprehensive literature research was conducted to search for modelling techniques that are easily applicable to phase change problems, as occur in the HyStEPs hybrid latent heat storage. The necessary theoretical framework to understand the development of the numerical model was outlined in Chapter 2.

Chapter 3 explained the assumptions and discretization made to model the phase change problem. The simple-to-use apparent heat capacity method (see Section 2.2.2) was chosen to account for latent heat implementation in the material. To discretize the energy equation (3.1), the finite element method was applied, using a rectangular mesh of four-noded bilinear elements to approximate the simulated domain. For the modelling of convection via the Navier Stokes equation (3.4), an existing finite di↵erence approximation method, documented in [36], was implemented, adapted to the problem and extended to fit the desired functionality. It should be noted though, that this code was used after unsuccessful attempts to discretize the Navier-Stokes equation by finite elements, using a penalty formulation for the pressure. In this case, the e↵ort to get a convergent numerical solution unfortunately failed. Nonetheless, this should not discourage fellow programmers from trying to apply this method for similar problems, as it was already done in [11].

A solid MATLAB implementation of the model explained in Chapter 3 was developed and documented in Chapter 4. This code can handle any desired layout of materials arranged on a rectangular domain, which is especially useful for PCM storage modelling with internal fins. The programming also aimed to achieve a high degree of vectorization to minimize the computational e↵ort. While simulations considering only heat conduction considered lead to very fast solutions, the sparsity patterns of the matrix equation systems to solve cannot be exploited the same way for simulations considering natural convection. For large numbers of elements, the large systems of equations to be solved decreases the computational performance significantly.

The simulation results could be successfully validated in three steps, as documented in Chapter 5. The modelling of heat conduction could be verified as very accurate, when compared to the analytical solution of the one-dimensional Stefan problem. Depending on the choice of computational parameters, relative errors of the evaluated quantities of well under 1 % could be obtained. Also, fast convergence of the solution with decreasing time step size and refinement of the mesh size were observed.

Regarding the validation of the convection model, a good agreement of the model predictions with the experimental results published in [5] was found, whilst the results of the mesh convergence analysis were not as promising as those for the one-dimensional heat conduction model.

Additionally, a cross-validation of the convection model, given by simulations run in the CFD software ANSYS Fluent, was presented in Section 5.3, which also showed only small discrepancies between the two simulation methods.

Finally, parameter studies conducted with the developed MATLAB model, are given in Chapter 6. Di↵erent cavity designs were analysed with respect to total enthalpy input per surface area on the steam storage, while sustaining a fast charging and discharging speed.

Di↵erent cavity orientations were studied with respect to their position on the steam storage, which, in the charging process, showed a significant influence of natural convection for upwards facing cavities and horizontal cavities, but diminishing influence for downwards facing cavities. During discharging, the e↵ect of natural convection was found to be insignificant, regardless of the cavity orientation.

## 7.2 Scope and limitations

The developed MATLAB model has proven to be an e↵ective tool for modelling phase change problems, considering both conduction and natural convection heat transfer, in typical PCM cavities. This was shown by conducting informative parameter studies, which provide an example of the methodology for further investigations.

The code can also be used as a basis for model order reduction, which will be carried out by fellow research colleagues. This will ultimately lead to a reduced order model, which is essential for the real-time operation and control of the HyStEPs hybrid latent heat storage.

In comparison to commercially available code, such as ANSYS Fluent, the developed MATLAB model provides high usability in regards to running and evaluating parameter studies, which can be run parallel on any desired number of sessions/computing units. Furthermore, modifications to the existing code can be easily made, aiding model order reduction by using MATLAB directly. This is a major benefit for future work carried out in the HyStEPs project.

Having said this, some issues remain with the current MATLAB model. The nature of the apparent heat capacity method brings some diculties in handling isothermal phase change problems. Although it could be obtained from Fluent simulations (see Section 6.2.1) that isothermal phase change can be approximated well by a very small mushy region, these simulations still require careful selection of computational parameters. As advised in the validation in Section 5.1, the magnitudes of the time step and mushy region size should ensure that the area of the mushy region progresses by overlapping its previous position. Otherwise, finite elements would skip the phase change and would not be exposed to the e↵ect of the latent heat. Even if these limitations are considered, an approximation error still remains, which is not the case for certain implementations of the enthalpy method, introduced in Section 2.2.1. Another option is source-based methods, in which the latent heat evolution is accounted for by the definition of a source term, instead of accounting for the latent heat in the specific heat coecient. Such an approach was followed for some time to further improve the model capabilities, though it could not completely be put into practice, for reasons of limited time during this work. However, it is highly suggested to take up this approach again, which is outlined in more detail in the next section (see Section 7.3.3).

Although the parameter studies reported in Chapter 6 were based on preliminary design and material choices of the HyStEPs cavity, at least the methodological approach can be used for the purpose of future parameter studies. The findings of Section 6.3, with regard to the e↵ect of natural convection in a typical PCM cavity, could prove to be a valuable insight for the development of the operating and control strategy in the HyStEPs project. For example, these results indicate a far less complicated charging behaviour for some cavity orientations, which implies a reduced number of sensors to be placed in the latter.

The observed e↵ects can of course occur very di↵erently for another PCM material or a modified cavity design. PCM materials with higher viscosity might lead to larger di↵erences between the simulated cavity orientations, as will be opposite for low viscosity substances. The aspect ratio of PCM cavities has a significant e↵ect on the heat transfer performances, as shown in reference [48], where, based on computational results, it was obtained that the melting rate increases as the aspect ratio of a rectangular cavity increases.

The reported e↵ect of natural convection during charging and the sucient approximation of the discharging process by only considering heat conduction was also reported in reference [27] for similar PCM cavity geometry.

## 7.3 Suggested future objectives

### 7.3.1 Subsequent parameter studies

As mentioned before, the design of the latent heat cavities, which will be mounted on a Ruths steam storage as part of the HyStEPs project, is not fixed at this time. Therefore, the parameter studies conducted in this thesis may not be particularly valuable for further investigations. Further parameter studies like the ones conducted in this thesis are suggested for the final design decision and to thoroughly analyse the characteristics of the chosen PCM cavity.

Building on the results of work, a more thorough parameter study is currently underway to provide a detailed foundation for design optimization of the considered PCM cell geometry [49]. Therein, varying aluminium proportions, varying fin spacings and di↵erent orientations of the PCM cell were simulated and their impact on charging/discharging speed was analysed.

### 7.3.2 Improvements of the existing MATLAB model

Some improvements to the developed MATLAB model may help reduce computational e↵ort while maintaining the same level of accuracy. One suggestion is to implement adaptive time stepping, like for example presented in [50], instead of a fixed time step. Choosing a lager time step automatically could help speed up simulations in which temperature gradients are small and fluid flow is slow.

Another option is to specify separate fixed time steps or adaptive time step schemes for the finite element heat conduction problem and the finite di↵erence Navier-Stokes solution method. Since the time step is generally limited by the CFL condition and not by the desired accuracy of the much faster FEM solution method, this would probably only lead to slightly better computational performance.

Furthermore, adaptive mesh or even moving mesh discretization methods, as documented in [11], could be considered for implementation. However, in this case, many of the faster preprocessing methods in the current code could note be used, and the benefit to performance is questionable

### 7.3.3 Fundamental modifications to the existing MATLAB model

#### Source-based method

As mentioned in Section 7.2, the implementation of a source-based method was attempted for some time during this work. Although the existing code was altered by the guidelines given in several references ([51] and [52]), the outcome was not satisfying since the simulation results were not reliable. For further reference, a brief explanation of the source-based method, is given here. This method could be especially important for accurately modelling isothermal phase change problems.

In the source-based method, the governing energy equation (2.5) for the apparent heat capacity method is written as

$$\frac{\partial(\rho cT)}{\partial t} = \nabla(k\nabla T) - \nabla(\rho cT \cdot \mathbf{u}) + S \ . \tag{7.1}$$

In this form however, *c* takes a fixed value instead of the apparent heat capacity, and the source term

$$S = -\delta H \frac{\delta a}{\delta t} \tag{7.2}$$

is added to the equation. Applying Galerkin's method to this equation, as done in Section 3.2.2, yields again a finite element system of the form

$$\begin{Bmatrix} C \end{Bmatrix} \begin{Bmatrix} \dot{T} \end{Bmatrix} + \begin{bmatrix} K \end{Bmatrix} \begin{Bmatrix} T \end{Bmatrix} = \begin{bmatrix} R \end{bmatrix} \begin{Bmatrix} \end{Bmatrix} \begin{Bmatrix} T \end{Bmatrix} \tag{3.23}$$

Here, the source term is incorporated in the system matrices [*K*] and [*R*], and the equation system (3.23) can be solved in each time step by an iterative scheme as highlighted in the following:


The obtained solution for the temperature field will not necessarily be consistent with the current liquid fraction field. For instance, a predicted nodal temperature *T* 6= *T<sup>m</sup>* would not be consistent with a liquid fraction *a* 2 (0*,* 1). That is why the value of the nodal liquid fraction field is updated so that on subsequent iterations of step 1 and 2 the predicted nodal temperature is 'driven' to *Tm*, until a desired accuracy is reached.

Besides the approximation of the source term, the key to the source-based iteration is the method of updating the liquid fraction *a*. Numerous schemes for this purpose are found in literature, of which the reference [51] is emphasized here.

#### Enthalpy method

A di↵erent approach to the phase change problem would be to consider enthalpy *H* as dependent variable, instead of the temperature *T*. Although this approach would require serious modifications to the existing MATLAB implementation, it is possible to build on this work and reuse the existing formalism. For further reference, a good example of a finite element implementation of the enthalpy method can be found in [9].

## Bibliography


## Appendix A

## MATLAB code listings

This appendix contains listings of the main functions of the developed MATLAB model.

Code Listing A.1: Code of main function PCMsolver

```
1 %% PCM Solver 2 % Author : Lukas Kasper
 3 % Last modified : 03.05.2019 4 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5
 6 function PCMsolver(input_data)
7 dbstop if error; 8 % clearvars -except input_data ; 9 % close all; 10
11 strTimestamp = datestr(now() ,30); tic; % get timestamp
12 fprintf( Starting simulation / timestamp % s \ n ,strTimestamp);
13
14 %% preprocessing
15 % load input variables 16 eval(input_data);
17 % additional input variables 18 non = param.mesh.non;
19 not = param.mesh.not;
20
21 % initialize waittime bar 22 if param.display.waitbar == 1
23 f = waitbar(0, Fortschritt ); timearray = 0; up = 0;
24 end 25
26 % mesh generation
27 [param] = mesh2d(param);
28 v2struct(param.geom);
29 v2struct(param.mesh);
30 v2struct(param.icbc);
31
32 % visualize material in domain 33 if param.display.matplot == 1
34 scatterbar3(center_coord(1,:) ,center_coord(2,:) ,... 35 elem_mat ,x_length/max(nx,ny));
36 end 37
38 % initialize arrays for data saving
39 Tt = zeros(non ,ceil(not/param.save.intervall));
40 s = 1;
41 if param.convection == 1
42 Ut = zeros(nx -1,ny,ceil(not/param.save.intervall));
43 Vt = zeros(nx ,ny -1,ceil(not/param.save.intervall));
44 end 45
46 % initialize nodal temperature matrix 47 T_new = zeros(non ,1);
48 T_old = T_init(:);
49 T_0 = zeros(not ,1);
50
51 % initialize staggered grid velocity matrices
52 U_old = zeros(nx -1,ny);
53 V_old = zeros(nx ,ny -1);
54 U_new = zeros(nx -1,ny);
55 V_new = zeros(nx ,ny -1);
56
57 % add boundaries to the velocity matriced
```

```
58 Uwb = [zeros(1,ny); U_old(:,:) ;zeros(1,ny)]; % east and west boundaries 59 Uwb = [-Uwb(:,1) Uwb -Uwb(:,end )]; % north and south boundaries 60 Vwb = [zeros(nx ,1) V_old (:,:) zeros(nx ,1)]; % north and south boundaries 61 Vwb = [-Vwb(1,:); Vwb ;-Vwb(end ,:)]; % east and west boundaries 62 Ue = avg(Uwb ) ; 63 Ve = avg(Vwb);
 64
 65 % initialize solid parts
 66 elempcm = elem_mat(reshape (1:nel ,nx,ny)) == 1;
 67 T_ij = reshape(T_old(:),nx+1,ny+1); % staggered grid temperatures
 68 %Tl = param .pcm.Tm+ param . pcm. dTm ; % temperature of complete
             liquidification
 69 Tl = param.pcm.Tm-param.pcm.dTm; % temperature of partial
             liquidification
 70 % p_el = phase of element [0 = solid , 1 = liquid ]
 71 p_el = false(nx ,ny);
 72
 73 % initialize heat flow and enthalpy arrays
 74 Q_in = zeros(not ,1);
 75 Q_out = zeros(not ,1);
 76 H = zeros(not ,1);
 77 H_l = zeros(not ,1);
 78
 79 %% FEM preprocessing
 80 % calculate shape matrices 81 [gauss] = shape_mat(param);
 82 % preprocess FEM data 83 [gauss] = FEM_preprocessor(T_old ,gauss ,param);
 84
 85 %% time stepping
 86 for t = 1:not 87
 88 %% FEM calculation 89 [C,K,R,Qin_e ,Qout_e ,H(t),H_l(t)] = ... 90 heat2Delem_new(t,T_old ,Uwb ,Vwb ,p_el ,gauss ,param);
 91
 92 % heat flows and enthalpy
 93 %%{
94 if t>1 95 Q_in(t) = Q_in(t-1) + Qin_e; % input heat flow 96 Q_out(t) = Q_out(t-1) + Qout_e; % output heat flow 97 else 98 Q_in(t) = Qin_e; % input heat flow 99 Q_out(t) = Qout_e; % output heat flow 100 end 101 % }
102
103 %% solve equation for temperatures at t+1 104 T_new = solve_euler_back(T_old ,C,K,R,param);
105
106 % refresh staggered grid temperatures
107 T_ij = reshape(T_new ,nx+1,ny+1); % staggered grid temperatures
108 if param.convection == 1
109 p_el = min(min( T_ij(1: end -1,1:ny),T_ij(1: end -1,2:ny+1)),... 110 min(T_ij(2: end ,1:ny),T_ij(2: end ,2:ny+1) ))>Tl... 111 & elempcm;
112 end 113
114 %% solve velocities at t+1 115 if any(any(p_el)) == 1
116
117 %% operating temperature for boussinesq term 118 % mean temperature in liquid region
119 T_0(t) = mean(sum(T_new(IEN(:,reshape(p_el ,[],1)==1)))/4);
120 % fixed temperature
121 % T_0 ( t) = 230;
122 param.pcm.T_0 = T_0(t);
123
124 %% solve navier stokes equations
125 [U_new ,V_new] = solve_navierstokes(U_old ,V_old ,T_ij ,p_el ,param);
126
127 %% reset velocities in solid part
128 U_new = U_new .*( p_el (1:nx -1,:) & p_el (2:nx ,:));
129 V_new = V_new .*( p_el (:,1:ny -1) & p_el (:,2:ny));
130
131 %% add boundary values to velocitiy matrices 132 % in domain phase boundary
133 U_new(:,1: end -1) = U_new(:,1: end -1)-U_new(:,2:end ).*... 134 (p_el(1: end -1,2: end)==1 & p_el(1: end -1,2: end )~=p_el(1: end -1,1: end -1));
```

```
135 U_new(:,2: end) = U_new(:,2: end)-U_new(:,1: end -1).*... 136 (p_el(1: end -1,1: end -1)==1 & p_el(1: end -1,1: end -1)~=p_el(1: end -1,2: end));
137 V_new(1: end -1,:) = V_new(1: end -1,:)-V_new(2: end ,:).*... 138 (p_el(2: end ,1: end -1)==1 & p_el(2: end ,1: end -1)~=p_el(1: end -1,1: end -1));
139 V_new(2: end ,:) = V_new(2: end ,:)-V_new(1: end -1,:).*... 140 (p_el(1: end -1,1: end -1)==1 & p_el(1: end -1,1: end -1)~=p_el(2: end ,1: end -1));
141
142 % domain boundaries 143 Uwb = [zeros(1,ny+2); [-U_new(:,1) U_new(:,:) -U_new(:, end )];... 144 zeros(1,ny+2)];
145 Vwb = [[0 -V_new (1,:) 0] ;[zeros(nx ,1) V_new (:,:) zeros(nx ,1)];... 146 [0 -V_new(end ,:) 0]];
147
148 else 149 U_new = U_old;
150 V_new = V_old;
151 end 152
153 %% waitbar update
154 if param.display.waitbar == 1
155 [f,timearray ,up] = waitbar2D_simple(f,param ,timearray ,toc ,t,up);
156 end 157 %% dump variables 158 if mod(t,param.save.intervall) == 0
159 Tt(:,s) = T_new;
160 if param.convection == 1
161 Ut(:,:,s) = U_new;
162 Vt(:,:,s) = V_new;
163 end 164 s = s+1;
165 end 166 if param.save.make_snapshot == 1
167 if mod(t,param.save.snapshot_int)==0
168 strFilename = sprintf( results / dump_ %s_t%03d . mat ,param.save.name
                          ,t);
169 if param.convection == 1
170 save(strFilename , C , K , R , t , T_new , T_old , U_new ,... 171 V_new , U_old , V_old , H , Q_in , Q_out , H_l , param ) ;
172 else 173 save(strFilename , C , K , R , t , T_new , T_old ,... 174 H , Q_in , Q_out , H_l , param ) ;
175 end 176 end 177 end 178
179 %% initialize T ,U , V for next timestep
180 T_old = T_new;
181 U_old = U_new;
182 V_old = V_new;
183 end 184
185 %% postprocessing
186 % delete waitbar 187 if param.display.waitbar == 1
188 delete(f);
189 end 190 % save data 191 strFilename = sprintf( results / dump_ %s_t%03 d_final .mat ,param.save.name ,t);
192 save(strFilename);
193
194 % plots after total time 195 if param.display.Tplot == 1
196 makeplot(T_new ,param);
197 end 198 if param.display.Vplot == 1
199 makeplot_konv(U_new ,V_new ,param);
200 end 201
202 % show simulation time 203 strTimestamp = datestr(now() ,30); % get timestamp
204 fprintf( Simulation complete / timestamp % s \n ,strTimestamp);
205 sim_time = toc;
206 fprintf( Total simulation time [ s]: %5.0 f \ n ,toc);
207 end
```
Code Listing A.2: Code of function mesh2d

```
1 %% creation of the 2 D mesh 2 function [param] = mesh2d(param)
 3 % The 2 D mesh is created . Coordinates and connectivities of nodes are 4 % computed . Boundary elements and element materials are flagged . 5 % All information is written to the structure array param . 6
 7 %% load input data 8 nx = param.mesh.nx;
 9 ny = param.mesh.ny;
10 non = param.mesh.non;
11 nel = param.mesh.nel;
12 v2struct(param.geom);
13
14 %% generate coordinates of nodes 15 x0 = linspace(0,x_length ,nx+1); % equal bisection of the x nodes 16 y0 = linspace(0,y_length ,ny+1); % equal bisection of the y nodes 17 x = zeros(1,non);
18 y = zeros(1,non);
19 for i = 1:(ny+1)
20 for j = 1:(nx+1)
21 x( (i-1)*(nx+1)+j ) = x0(j);
22 y( (i-1)*(nx+1)+j ) = y0(i);
23 end 24 end 25 param.mesh.x = x;
26 param.mesh.y = y;
27
28 %% generate the IEN connectivity array
29 % IEN = identify element nodes 30 IEN = zeros(4,nel);
31 rowcount = 0;
32 for elementcount = 1:nel 33 IEN(1,elementcount) = elementcount + rowcount;
34 IEN(2,elementcount) = elementcount + 1 + rowcount;
35 IEN(3,elementcount) = elementcount + (nx + 2) + rowcount;
36 IEN(4,elementcount) = elementcount + (nx + 1) + rowcount;
37 if mod(elementcount ,nx) == 0
38 rowcount = rowcount + 1;
39 end 40 end 41 param.mesh.IEN = IEN;
42
43 %% flag boundary elements 44 flag = zeros(nel ,2);
45 for e = 1:nel 46 if mod(e,nx) == 1
47 flag(e,1) = 3; % 3 = left side element 48 elseif mod(e,nx) == 0
49 flag(e,1) = 1; % 1 = right side element 50 end 51 for i = 1:nx 52 flag(i,2) = 4; % 4 = lower side element 53 flag((nel -i+1) ,2) = 2; % 2 = upper side element 54 end 55 end 56 param.mesh.flag = flag;
57
58 %% flag element matieral 59 % coordinates of element center 60 center_coord = zeros(2,nel);
61 center_coord(:,:) = [(x(IEN(2,:))+x(IEN(1,:)))/2;(y(IEN(3,:))+y(IEN(2,:)))
            /2];
62
63 % flag element matieral 64 elem_mat = zeros(nel ,1);
65 for i = length(x_length_vec):-1:1
66 for j = length(y_length_vec):-1:1
67 x = sum(x_length_vec(1:i));
68 y = sum(y_length_vec(1:j));
69 elem_mat(center_coord(1,:)<x & center_coord(2,:)<y) = mat_mtrx(j,i);
70 end 71 end 72 param.mesh.elem_mat = elem_mat;
73 param.mesh.center_coord = center_coord;
74
75 %% generate staggered grid coordinates for convection 76 stag = zeros(param.mesh.nel ,2);
77 for e = 1:param.mesh.nel
```

```
78 if mod(e,param.mesh.nx) == 0
79 stag(e,1) = param.mesh.nx;
80 else 81 stag(e,1) = mod(e,param.mesh.nx);
82 end 83 stag(e,2) = ceil(e/param.mesh.nx);
84 end 85 param.mesh.stag = stag;
```
Code Listing A.3: Code of function shape mat

```
1 %% shape function calculation 2 function [gauss] = shape_mat(param)
 3 % If all elements are the same size , use this function to calculate the 4 % shape functions , derivate of shape functions and Jacobian to be later 5 % used in the FEM Calculation heat2Delem_2 . 6
 7 nodes = param.mesh.IEN(:,1);
 8 coord = [param.mesh.x(nodes); param.mesh.y(nodes)] ; 9 gp = [-0.57735027, 0.57735027 ]; % gauss points (2 Point Gauss - Legendre )
10
11 % for element terms 12 gauss.N1 = NmatHeat2D(gp(1),gp(1));
13 gauss.N1_q = gauss.N1 * gauss.N1;
14 [gauss.B1, gauss.detJ] = BmatHeat2D(gp(1),gp(1),coord);
15 gauss.B1_q = gauss.B1 * gauss.B1;
16
17 gauss.N2 = NmatHeat2D(gp(2),gp(1));
18 gauss.N2_q = gauss.N2 * gauss.N2;
19 gauss.B2 = BmatHeat2D(gp(2),gp(1),coord);
20 gauss.B2_q = gauss.B2 * gauss.B2;
21
22 gauss.N3 = NmatHeat2D(gp(2),gp(2));
23 gauss.N3_q = gauss.N3 * gauss.N3;
24 gauss.B3 = BmatHeat2D(gp(2),gp(2),coord);
25 gauss.B3_q = gauss.B3 * gauss.B3;
26
27 gauss.N4 = NmatHeat2D(gp(1),gp(2));
28 gauss.N4_q = gauss.N4 * gauss.N4;
29 gauss.B4 = BmatHeat2D(gp(1),gp(2),coord);
30 gauss.B4_q = gauss.B4 * gauss.B4;
31
32 % for boundary terms 33 eta = 1;
34 gauss.Ne1 = NmatHeat2D(eta ,gp(1));
35 gauss.Ne1_q = gauss.Ne1 * gauss.Ne1;
36 gauss.Ne2 = NmatHeat2D(eta ,gp(2));
37 gauss.Ne2_q = gauss.Ne2 * gauss.Ne2;
38 gauss.dS_e = 0.25*[eta -1 -eta -1 1+eta 1-eta]*coord(:,2);
39 eta = -1;
40 gauss.Nw1 = NmatHeat2D(eta ,gp(1));
41 gauss.Nw1_q = gauss.Nw1 * gauss.Nw1;
42 gauss.Nw2 = NmatHeat2D(eta ,gp(2));
43 gauss.Nw2_q = gauss.Nw2 * gauss.Nw2;
44 gauss.dS_w = - 0.25*[eta -1 -eta -1 1+eta 1-eta]*coord(:,2);
45
46 % for convection terms 47 gp = [-1,0,1]; % convection gauss points
48 w = [1/3,3/4,1/3]; % convection gauss weights
49 gauss.Nc1 = NmatHeat2D(gp(1),gp(1)); % eta -1, psi -1 50 [B1, detJ] = BmatHeat2D(gp(1),gp(1),coord);
51 mat_const = param.pcm.rho*param.pcm.c_l*detJ;
52 gauss.NBc1u = w(1)*w(1)*mat_const*gauss.Nc1 * B1(1,:);
53 gauss.NBc1v = w(1)*w(1)*mat_const*gauss.Nc1 * B1(2,:);
54
55 gauss.Nc2 = NmatHeat2D(gp(1),gp(2)); % eta -1, psi 0 56 B2 = BmatHeat2D(gp(1),gp(2),coord);
57 gauss.NBc2u = w(1)*w(2)*mat_const*gauss.Nc2 * B2(1,:);
58 gauss.NBc2v = w(1)*w(2)*mat_const*gauss.Nc2 * B2(2,:);
59
60 gauss.Nc3 = NmatHeat2D(gp(1),gp(3)); % eta -1, psi 1 61 B3 = BmatHeat2D(gp(1),gp(3),coord);
62 gauss.NBc3u = w(1)*w(3)*mat_const*gauss.Nc3 * B3(1,:);
63 gauss.NBc3v = w(1)*w(3)*mat_const*gauss.Nc3 * B3(2,:);
64
65 gauss.Nc4 = NmatHeat2D(gp(2),gp(1)); % eta 0 , psi -1 66 B4 = BmatHeat2D(gp(2),gp(1),coord);
```

```
67 gauss.NBc4u = w(2)*w(1)*mat_const*gauss.Nc4 * B4(1,:);
68 gauss.NBc4v = w(2)*w(1)*mat_const*gauss.Nc4 * B4(2,:);
69
70 gauss.Nc5 = NmatHeat2D(gp(2),gp(2)); % eta 0 , psi 0 71 B5 = BmatHeat2D(gp(2),gp(2),coord);
72 gauss.NBc5u = w(2)*w(2)*mat_const*gauss.Nc5 * B5(1,:);
73 gauss.NBc5v = w(2)*w(2)*mat_const*gauss.Nc5 * B5(2,:);
74
75 gauss.Nc6 = NmatHeat2D(gp(2),gp(3)); % eta 0 , psi 1 76 B6 = BmatHeat2D(gp(2),gp(3),coord);
77 gauss.NBc6u = w(2)*w(3)*mat_const*gauss.Nc6 * B6(1,:);
78 gauss.NBc6v = w(2)*w(3)*mat_const*gauss.Nc6 * B6(2,:);
79
80 gauss.Nc7 = NmatHeat2D(gp(3),gp(1)); % eta 1 , psi -1 81 B7 = BmatHeat2D(gp(3),gp(1),coord);
82 gauss.NBc7u = w(3)*w(1)*mat_const*gauss.Nc7 * B7(1,:);
83 gauss.NBc7v = w(3)*w(1)*mat_const*gauss.Nc7 * B7(2,:);
84
85 gauss.Nc8 = NmatHeat2D(gp(3),gp(2)); % eta 1 , psi 0 86 B8 = BmatHeat2D(gp(3),gp(2),coord);
87 gauss.NBc8u = w(3)*w(2)*mat_const*gauss.Nc8 * B8(1,:);
88 gauss.NBc8v = w(3)*w(2)*mat_const*gauss.Nc8 * B8(2,:);
89
90 gauss.Nc9 = NmatHeat2D(gp(3),gp(3)); % eta 1 , psi 1 91 B9 = BmatHeat2D(gp(3),gp(3),coord);
92 gauss.NBc9u = w(3)*w(3)*mat_const*gauss.Nc9 * B9(1,:);
93 gauss.NBc9v = w(3)*w(3)*mat_const*gauss.Nc9 * B9(2,:);
94
95 end
```
Code Listing A.4: Code of shape matrix function NmatHeat2D

```
1 %% N matrix Heat2D 2 function N = NmatHeat2D(eta ,psi)
3 % Computes shape functions matrix N 4
5 N = 0.25 * [(1-psi)*(1-eta) (1-psi)*(1+eta) ... 6 (1+psi)*(1+eta) (1+psi)*(1-eta)]; % shape functions
```
Code Listing A.5: Code of shape matrix function BmatHeat2D

```
1 %% B matrix Heat2D 2 function [B, detJ] = BmatHeat2D(eta ,psi ,coord)
 3 % Computes B matrix 4 % = derivative of the shape functions matrix N 5
 6 % Calculate the Grad ( N ) matrix 7 GN = 0.25 * [psi -1 1-psi 1+psi -psi -1;
 8 eta -1 -eta -1 1+eta 1-eta];
 9
10 J = GN*coord; % compute Jacobian matrix 11 detJ = det(J); % Jacobian 12
13 B = J\GN; % compute the B matrix
```
Code Listing A.6: Code of FEM preprocessing function FEM preprocessor

```
1 %% FEM preprocessor
2 function [gauss] = FEM_preprocessor(T_old ,gauss ,param)
3 % preprocesses element properties for FEM calculation 4
5 % load important variables 6 elem_mat = param.mesh.elem_mat;
7 IEN = param.mesh.IEN;
8 nel = param.mesh.nel;
9 elem_pcm = find(elem_mat==1); % element indices of PCM 10
11 %% intialize preprocessed data arrays
12 c1 = zeros(nel ,1);
13 c2=c1;c3=c1;c4=c1;k1=c1;k2=c1;k3=c1;k4=c1;rho=c1;
14 flag_el=c1;calc_h=c1;
```

```
15
16 %% preprocess data for 4 gauss points per element 17
18 % heat capacities and heat conductivities 19 % gauss point 1 20 T1 = gauss.N1*T_old(IEN(:,:));
21 c1(elem_mat ==1) = c_var(T1(elem_mat ==1),param.pcm);
22 k1(elem_mat ==1) = k_var(T1(elem_mat ==1),param.pcm);
23 for m = 2:size(param.geom.material ,1)
24 str1 = strcat( param . ,param.geom.material{m,2}, .c_s ) ;
25 str2 = strcat( param . ,param.geom.material{m,2}, .k_s ) ;
26 c1(elem_mat==m) = eval(str1);
27 k1(elem_mat==m) = eval(str2);
28 end 29 % gauss point 2 30 T2 = gauss.N2*T_old(IEN(:,:));
31 c2(elem_mat ==1) = c_var(T2(elem_mat ==1),param.pcm);
32 k2(elem_mat ==1) = k_var(T2(elem_mat ==1),param.pcm);
33 for m = 2:size(param.geom.material ,1)
34 str1 = strcat( param . ,param.geom.material{m,2}, .c_s ) ;
35 str2 = strcat( param . ,param.geom.material{m,2}, .k_s ) ;
36 c2(elem_mat==m) = eval(str1);
37 k2(elem_mat==m) = eval(str2);
38 end 39 % gauss point 3 40 T3 = gauss.N3*T_old(IEN(:,:));
41 c3(elem_mat ==1) = c_var(T3(elem_mat ==1),param.pcm);
42 k3(elem_mat ==1) = k_var(T3(elem_mat ==1),param.pcm);
43 for m = 2:size(param.geom.material ,1)
44 str1 = strcat( param . ,param.geom.material{m,2}, .c_s ) ;
45 str2 = strcat( param . ,param.geom.material{m,2}, .k_s ) ;
46 c3(elem_mat==m) = eval(str1);
47 k3(elem_mat==m) = eval(str2);
48 end 49 % gauss point 4 50 T4 = gauss.N4*T_old(IEN(:,:));
51 c4(elem_mat ==1) = c_var(T4(elem_mat ==1),param.pcm);
52 k4(elem_mat ==1) = k_var(T4(elem_mat ==1),param.pcm);
53 for m = 2:size(param.geom.material ,1)
54 str1 = strcat( param . ,param.geom.material{m,2}, .c_s ) ;
55 str2 = strcat( param . ,param.geom.material{m,2}, .k_s ) ;
56 c4(elem_mat==m) = eval(str1);
57 k4(elem_mat==m) = eval(str2);
58 end 59
60 % density of element 61 for m = 1:size(param.geom.material ,1)
62 str = strcat( param . ,param.geom.material{m,2}, .rho ) ;
63 rho(elem_mat==m) = eval(str);
64 end 65
66 % constants for enthalpy calculation 67 for m = 1:size(param.geom.material ,1)
68 str1 = strcat( param . ,param.geom.material{m,2}, .rho ) ;
69 str2 = strcat( param . ,param.geom.material{m,2}, .c_s ) ;
70 calc_h(elem_mat==m) = eval(str1)*gauss.detJ*eval(str2);
71 end 72
73 %% preprocess boundary information 74 % flag boundary elements by one identifier 75 flag_el(param.mesh.flag(:,1) == 3) = 3; % left side element 76 flag_el(param.mesh.flag(:,1) == 1) = 1; % right side element 77 flag_el(param.mesh.flag(:,2) == 2) = 2; % upper side element 78 flag_el(param.mesh.flag(:,2) == 4) = 4; % lower side element 79 flag_el((param.mesh.flag(:,1) == 3)&(param.mesh.flag(:,2) == 4) ) = 34;
80 flag_el((param.mesh.flag(:,1) == 3)&(param.mesh.flag(:,2) == 2) ) = 32;
81 flag_el((param.mesh.flag(:,1) == 1)&(param.mesh.flag(:,2) == 4) ) = 14;
82 flag_el((param.mesh.flag(:,1) == 1)&(param.mesh.flag(:,2) == 2) ) = 12;
83
84 gauss.west = find(param.mesh.flag(:,1) == 3);
85 gauss.east = find(param.mesh.flag(:,1) == 1);
86 yr=param.mesh.y(param.mesh.IEN(:,:));
87 yr1 = yr(1,:)+(yr(4,:)-yr(1,:))*(1 -0.57735027)/2;
88 yr2 = yr(1,:)+(yr(4,:)-yr(1,:))*(1+0.57735027) /2;
89 gauss.Tout1 = @(t) param.icbc.Tout(t,yr1);
90 gauss.Tout2 = @(t) param.icbc.Tout(t,yr2);
91 gauss.alpha_out1 = @(t) param.icbc.alpha_out(t,yr1);
92 gauss.alpha_out2 = @(t) param.icbc.alpha_out(t,yr2);
93
94 yl=param.mesh.y(param.mesh.IEN(:,:));
```

```
95 yl1 = yl(1,:)+(yl(4,:)-yl(1,:))*(1 -0.57735027)/2;
 96 yl2 = yl(1,:)+(yl(4,:)-yl(1,:))*(1+0.57735027) /2;
 97 gauss.Tin1 = @(t) param.icbc.Tin(t,yl1);
 98 gauss.Tin2 = @(t) param.icbc.Tin(t,yl2);
 99 gauss.alpha_in1 = @(t) param.icbc.alpha_in(t,yl1);
100 gauss.alpha_in2 = @(t) param.icbc.alpha_in(t,yl2);
101
102 %% preprocess sparse assembling relation 103 sp_mat_e = cell(nel ,1);
104 sp_vec_e = cell(nel ,1);
105 for e = 1:nel 106 je = param.mesh.IEN(:,e);
107
108 k = [je;je;je;je];
109 l = [je(1);je(1);je(1);je(1);... 110 je(2);je(2);je(2);je(2);... 111 je(3);je(3);je(3);je(3);... 112 je(4);je(4);je(4);je(4)];
113
114 sp_mat_e{e} = [k,l];
115 sp_vec_e{e} = [je ,[1;1;1;1]];
116 end 117 sp_mat = cell2mat(sp_mat_e);
118 sp_vec = cell2mat(sp_vec_e);
119
120 %% load variables to struct variable gauss
121 gauss.rho = rho;
122 gauss.elem_pcm = elem_pcm;
123 gauss.sp_mat = sp_mat;
124 gauss.sp_vec = sp_vec;
125 gauss.flag_el = flag_el;
126 gauss.calc_h = calc_h;
127 gauss.c1 = c1;
128 gauss.c2 = c2;
129 gauss.c3 = c3;
130 gauss.c4 = c4;
131 gauss.k1 = k1;
132 gauss.k2 = k2;
133 gauss.k3 = k3;
134 gauss.k4 = k4;
135
136 end
```
Code Listing A.7: Code of FEM matrix assembly function heat2Delem new

```
1 %% heat 2 D element new 2 function [C,K,R,Qin ,Qout ,H,H_l] = heat2Delem_new(t,T,U,V,p_el ,gauss ,param)
3 % Computes the finite element matrices 4
5 %% initialize variables and matrix cells 6 nel = param.mesh.nel;
7 non = param.mesh.non;
8 nx = param.mesh.nx;
9 elem_pcm = gauss.elem_pcm;
10 not_pcm = param.mesh.elem_mat ~= 1;
11 sp_mat = gauss.sp_mat;
12 sp_vec = gauss.sp_vec;
13 flag_el = gauss.flag_el;
14 calc_h = gauss.calc_h;
15 calc_h1 = gauss.calc_h;
16 calc_h2 = gauss.calc_h;
17 calc_h3 = gauss.calc_h;
18 calc_h4 = gauss.calc_h;
19 hl_1 = gauss.calc_h;
20 hl_2 = gauss.calc_h;
21 hl_3 = gauss.calc_h;
22 hl_4 = gauss.calc_h;
23 c1 = gauss.c1;
24 c2 = gauss.c2;
25 c3 = gauss.c3;
26 c4 = gauss.c4;
27 k1 = gauss.k1;
28 k2 = gauss.k2;
29 k3 = gauss.k3;
30 k4 = gauss.k4;
31
32 p_el_FEM = reshape(p_el ,[],1);
33 %C_e = cell (nel ,1);
```

```
34 %K_e = cell (nel ,1) ; 35 %R_e = cell (nel ,1) ; 36
 37 %% vectorized preprocessing of element data 38 % temperatures at gauss points
 39 T1 = gauss.N1*T(param.mesh.IEN(:,:));
 40 T2 = gauss.N2*T(param.mesh.IEN(:,:));
 41 T3 = gauss.N3*T(param.mesh.IEN(:,:));
 42 T4 = gauss.N4*T(param.mesh.IEN(:,:));
 43 Te1 = gauss.Ne1*T(param.mesh.IEN(:,:));
 44 Te2 = gauss.Ne2*T(param.mesh.IEN(:,:));
 45 Tw1 = gauss.Nw1*T(param.mesh.IEN(:,:));
 46 Tw2 = gauss.Nw2*T(param.mesh.IEN(:,:));
 47
 48 % for conductive terms 49 % heat capacities and heat conductivities at gauss points of pcm elements 50 c1(elem_pcm) = c_var(T1(elem_pcm),param.pcm);
 51 k1(elem_pcm) = k_var(T1(elem_pcm),param.pcm);
 52 c2(elem_pcm) = c_var(T2(elem_pcm),param.pcm);
 53 k2(elem_pcm) = k_var(T2(elem_pcm),param.pcm);
 54 c3(elem_pcm) = c_var(T3(elem_pcm),param.pcm);
 55 k3(elem_pcm) = k_var(T3(elem_pcm),param.pcm);
 56 c4(elem_pcm) = c_var(T4(elem_pcm),param.pcm);
 57 k4(elem_pcm) = k_var(T4(elem_pcm),param.pcm);
 58 % constants for enthalpy calculation at gauss points of pcm elements 59 [calc_h1(elem_pcm),hl_1(elem_pcm)] = h_var(T1(elem_pcm),param ,gauss.detJ);
 60 [calc_h2(elem_pcm),hl_2(elem_pcm)] = h_var(T2(elem_pcm),param ,gauss.detJ);
 61 [calc_h3(elem_pcm),hl_3(elem_pcm)] = h_var(T3(elem_pcm),param ,gauss.detJ);
 62 [calc_h4(elem_pcm),hl_4(elem_pcm)] = h_var(T4(elem_pcm),param ,gauss.detJ);
 63
 64 % for convective terms 65 % velocities at gauss points of pcm elements 66 if any(p_el_FEM) == 1
 67 ie = param.mesh.stag(:,1);
 68 je = param.mesh.stag(:,2);
 69 % eta -1, psi -1 70 u1 = (U(sub2ind(size(U),ie,je))+U(sub2ind(size(U),ie,je+1)))/2;
 71 v1 = (V(sub2ind(size(V),ie,je))+V(sub2ind(size(V),ie+1,je)))/2;
 72 % eta -1, psi 0 73 u2 = U(sub2ind(size(U),ie,je+1));
 74 v2 = (V(sub2ind(size(V),ie,je))+V(sub2ind(size(V),ie,je+1))+... 75 V(sub2ind(size(V),ie+1,je))+V(sub2ind(size(V),ie+1,je+1)) )/2;
 76 % eta -1, psi 1 77 u3 = (U(sub2ind(size(U),ie,je+1))+U(sub2ind(size(U),ie,je+2)))/2;
 78 v3 = (V(sub2ind(size(V),ie,je+1))+V(sub2ind(size(V),ie+1,je+1)))/2;
 79 % eta 0 , psi -1 80 u4 = (U(sub2ind(size(U),ie,je))+U(sub2ind(size(U),ie,je+1))+... 81 U(sub2ind(size(U),ie+1,je))+U(sub2ind(size(U),ie+1,je+1)) )/2;
 82 v4 = V(sub2ind(size(V),ie,je+1));
 83 % eta 0 , psi 0 84 u5 = (U(sub2ind(size(U),ie,je+1))+U(sub2ind(size(U),ie+1,je+1)))/2;
 85 v5 = (V(sub2ind(size(V),ie+1,je))+V(sub2ind(size(V),ie+1,je+1)))/2;
 86 % eta 0 , psi 1 87 u6 = (U(sub2ind(size(U),ie,je+1))+U(sub2ind(size(U),ie,je+2))+... 88 U(sub2ind(size(U),ie+1,je+1))+U(sub2ind(size(U),ie+1,je+2)) )/2;
 89 v6 = V(sub2ind(size(V),ie+1,je+1));
 90 % eta 1 , psi -1 91 u7 = (U(sub2ind(size(U),ie+1,je))+U(sub2ind(size(U),ie+1,je+1)))/2;
 92 v7 = (V(sub2ind(size(V),ie+1,je))+V(sub2ind(size(V),ie+2,je)))/2;
 93 % eta 1 , psi 0 94 u8 = U(sub2ind(size(U),ie+1,je+1));
 95 v8 = (V(sub2ind(size(V),ie+1,je))+V(sub2ind(size(V),ie+1,je+1))+... 96 V(sub2ind(size(V),ie+2,je))+V(sub2ind(size(V),ie+2,je+1)) )/2;
 97 % eta 1 , psi 1 98 u9 = (U(sub2ind(size(U),ie+1,je+1))+U(sub2ind(size(U),ie+1,je+2)))/2;
 99 v9 = (V(sub2ind(size(V),ie+1,je+1))+V(sub2ind(size(V),ie+2,je+1)))/2;
100
101 %% velocity correction 102 % mark elements for velocity correction 103 correct_u = zeros(nel ,1);
104 correct_v = zeros(nel ,1);
105 solid=p_el_FEM==0;
106
107 % middle element 108 correct_u(flag_el==0) = solid(circshift(flag_el==0,1))*1+... 109 solid(circshift(flag_el==0,-1))*(-1)+... 110 ( solid(circshift(flag_el==0,1)) & ... 111 solid(circshift(flag_el==0,-1)))*2;
112 correct_v(flag_el ==0) = solid(circshift(flag_el==0,nx))*1+... 113 solid(circshift(flag_el==0,-nx))*(-1)+...
```
 ( solid(circshift(flag\_el==0,nx)) & ... <sup>115</sup> solid(circshift(flag\_el==0,-nx)) )\*2; % left side element <sup>117</sup> correct\_u(flag\_el ==3) = -1+solid(circshift(flag\_el ==3,1))\*3; correct\_v(flag\_el ==3) = solid(circshift(flag\_el==3,nx))\*1+... <sup>119</sup> solid(circshift(flag\_el==3,-nx))\*(-1)+... <sup>120</sup> ( solid(circshift(flag\_el==3,nx)) & ... <sup>121</sup> solid(circshift(flag\_el==3,-nx)) )\*2; correct\_u(flag\_el ==32) = -1+solid(circshift(flag\_el ==32,1))\*3; correct\_v(flag\_el ==32) = 1+solid(circshift(flag\_el==32,-nx))\*1; correct\_u(flag\_el ==34) = -1+solid(circshift(flag\_el ==34,1))\*3; correct\_v(flag\_el ==34) = -1+solid(circshift(flag\_el==32,nx))\*3; % right side element <sup>127</sup> correct\_u(flag\_el ==1) = 1+solid(circshift(flag\_el==1,-1))\*1; correct\_v(flag\_el ==1) = solid(circshift(flag\_el==1,nx))\*1+... <sup>129</sup> solid(circshift(flag\_el==1,-nx))\*(-1)+... <sup>130</sup> ( solid(circshift(flag\_el==1,nx)) & ... <sup>131</sup> solid(circshift(flag\_el==1,-nx)) )\*2; correct\_u(flag\_el ==12) = 1+solid(circshift(flag\_el==12,-1))\*1; correct\_v(flag\_el ==12) = 1+solid(circshift(flag\_el==12,-nx))\*1; correct\_u(flag\_el ==14) = 1+solid(circshift(flag\_el==14,-1))\*1; correct\_v(flag\_el ==14) = -1+solid(circshift(flag\_el==14,nx))\*3; % upper side element <sup>137</sup> correct\_u(flag\_el==2) = solid(circshift(flag\_el==2,1))\*1+... <sup>138</sup> solid(circshift(flag\_el==2,-1))\*(-1)+... <sup>139</sup> ( solid(circshift(flag\_el==2,1)) & ... <sup>140</sup> solid(circshift(flag\_el==2,-1)) )\*2; correct\_v(flag\_el ==2) = 1+solid(circshift(flag\_el==2,-nx))\*1; % lower side element <sup>143</sup> correct\_u(flag\_el==4) = solid(circshift(flag\_el==4,1))\*1+... <sup>144</sup> solid(circshift(flag\_el==4,-1))\*(-1)+... <sup>145</sup> ( solid(circshift(flag\_el==4,1)) & ... <sup>146</sup> solid(circshift(flag\_el==4,-1)) )\*2; correct\_v(flag\_el ==4) = -1+solid(circshift(flag\_el==4,nx))\*3; % perform velocity correction <sup>150</sup> u=[u1 ,u2 ,u3 ,u4 ,u5 ,u6 ,u7 ,u8 ,u9]; v=[v1 ,v2 ,v3 ,v4 ,v5 ,v6 ,v7 ,v8 ,v9]; u(correct\_u==2,:) = 0; u(correct\_u ==1,:) = max(0,u(correct\_u ==1,:)); u(correct\_u==-1,:) = min(0,u(correct\_u==-1,:)); v(correct\_v==2,:) = 0; v(correct\_v ==1,:) = max(0,v(correct\_v ==1,:)); v(correct\_v==-1,:) = min(0,v(correct\_v==-1,:)); u1 = u(:,1); u2 = u(:,2); u3 = u(:,3); u4 = u(:,4); u5 = u(:,5); u6 = u(:,6); u7 = u(:,7); u8 = u(:,8); u9 = u(:,9); v1 = v(:,1); v2 = v(:,2); v3 = v(:,3); v4 = v(:,4); v5 = v(:,5); v6 = v(:,6); v7 = v(:,7); v8 = v(:,8); v9 = v(:,9); end <sup>163</sup> %% calculate element matrices <sup>165</sup> Ce2 = gauss.detJ\*(... <sup>166</sup> bsxfun(@times ,gauss.N1\_q ,permute(gauss.rho .\* c1 ,[3 1 2]))+... <sup>167</sup> bsxfun(@times ,gauss.N2\_q ,permute(gauss.rho .\* c2 ,[3 1 2]))+... <sup>168</sup> bsxfun(@times ,gauss.N3\_q ,permute(gauss.rho .\* c3 ,[3 1 2]))+... <sup>169</sup> bsxfun(@times ,gauss.N4\_q ,permute(gauss.rho .\* c4 ,[3 1 2]))); C\_e = reshape(Ce2 ,4\*4\*nel ,1); KCe = gauss.detJ\*(... <sup>173</sup> bsxfun (@times ,gauss.B1\_q ,permute(k1 ,[3 1 2]))+... <sup>174</sup> bsxfun (@times ,gauss.B2\_q ,permute(k2 ,[3 1 2]))+... <sup>175</sup> bsxfun (@times ,gauss.B3\_q ,permute(k3 ,[3 1 2]))+... <sup>176</sup> bsxfun (@times ,gauss.B4\_q ,permute(k4 ,[3 1 2]))); right = param.mesh.flag(:,1) == 1; % right side of PCM <sup>179</sup> KB1e = gauss.dS\_e\*( ... <sup>180</sup> bsxfun(@times ,gauss.Ne1\_q ,permute(right .\* gauss.alpha\_out1(t) ,[3 1 2]))+... bsxfun(@times ,gauss.Ne2\_q ,permute(right .\* gauss.alpha\_out2(t) ,[3 1 2]))); RB1e = gauss.dS\_e\*( ... <sup>184</sup> bsxfun(@times ,gauss.Ne1 , permute(right .\* gauss.alpha\_out1(t).\*gauss. Tout1(t) ,[3 1 2]))+... bsxfun(@times ,gauss.Ne2 , permute(right .\* gauss.alpha\_out2(t).\*gauss. Tout2(t) ,[3 1 2]))); left = param.mesh.flag(:,1) == 3; % left side of PCM <sup>188</sup> KB3e = -gauss.dS\_w\*( ...

```
189 bsxfun(@times ,gauss.Nw1_q ,permute(left .* gauss.alpha_in1(t) ,[3 1 2])) +... 190 bsxfun(@times ,gauss.Nw2_q ,permute(left .* gauss.alpha_in2(t) ,[3 1 2]))
                                                   ) ;
191
192 RB3e = -gauss.dS_w*( ... 193 bsxfun(@times ,gauss.Nw1 , permute(left .* gauss.alpha_in1(t).*gauss.
Tin1(t) ,[3 1 2]))+... 194 bsxfun(@times ,gauss.Nw2 , permute(left .* gauss.alpha_in2(t).*gauss.
                                                    Tin2(t) ,[3 1 2])));
195
196 % convection term 197 if any(p_el_FEM) == 1
198 KCe = KCe + bsxfun(@times ,gauss.NBc1u ,permute(p_el_FEM .* u1 ,[3 1 2]))+ ... 199 bsxfun(@times ,gauss.NBc2u ,permute(p_el_FEM .* u2 ,[3 1 2]))+ ... 200 bsxfun(@times ,gauss.NBc3u ,permute(p_el_FEM .* u3 ,[3 1 2]))+ ... 201 bsxfun(@times ,gauss.NBc4u ,permute(p_el_FEM .* u4 ,[3 1 2]))+ ... 202 bsxfun(@times ,gauss.NBc5u ,permute(p_el_FEM .* u5 ,[3 1 2]))+ ... 203 bsxfun(@times ,gauss.NBc6u ,permute(p_el_FEM .* u6 ,[3 1 2]))+ ... 204 bsxfun(@times ,gauss.NBc7u ,permute(p_el_FEM .* u7 ,[3 1 2]))+ ... 205 bsxfun(@times ,gauss.NBc8u ,permute(p_el_FEM .* u8 ,[3 1 2]))+ ... 206 bsxfun(@times ,gauss.NBc9u ,permute(p_el_FEM .* u9 ,[3 1 2]))+ ... 207 bsxfun(@times ,gauss.NBc1v ,permute(p_el_FEM .* v1 ,[3 1 2]))+ ... 208 bsxfun(@times ,gauss.NBc2v ,permute(p_el_FEM .* v2 ,[3 1 2]))+ ... 209 bsxfun(@times ,gauss.NBc3v ,permute(p_el_FEM .* v3 ,[3 1 2]))+ ... 210 bsxfun(@times ,gauss.NBc4v ,permute(p_el_FEM .* v4 ,[3 1 2]))+ ... 211 bsxfun(@times ,gauss.NBc5v ,permute(p_el_FEM .* v5 ,[3 1 2]))+ ... 212 bsxfun(@times ,gauss.NBc6v ,permute(p_el_FEM .* v6 ,[3 1 2]))+ ... 213 bsxfun(@times ,gauss.NBc7v ,permute(p_el_FEM .* v7 ,[3 1 2]))+ ... 214 bsxfun(@times ,gauss.NBc8v ,permute(p_el_FEM .* v8 ,[3 1 2]))+ ... 215 bsxfun(@times ,gauss.NBc9v ,permute(p_el_FEM .* v9 ,[3 1 2]));
216 end 217 K_e = reshape(KCe + KB1e + KB3e ,4*4*nel ,1);
218 R_e = reshape(RB1e + RB3e ,4*nel ,1);
219
220 %% assemble to global matrices and vectors 221 C = sparse(sp_mat(:,1),sp_mat(:,2),C_e ,non ,non);
222 K = sparse(sp_mat(:,1),sp_mat(:,2),K_e ,non ,non);
223 R = sparse(sp_vec(:,1),sp_vec(:,2),R_e ,non ,1);
224
225 %% calculate enthalpy values 226 H = sum( calc_h1(elem_pcm)+calc_h2(elem_pcm)+... 227 calc_h3(elem_pcm)+calc_h4(elem_pcm)) + ... 228 sum( calc_h(not_pcm) .*(T1(not_pcm)+... 229 T2(not_pcm)+T3(not_pcm)+T4(not_pcm)));
230
231 H_l = sum( hl_1(elem_pcm)+hl_2(elem_pcm)+hl_3(elem_pcm)+hl_4(elem_pcm) );
232
233 Qout = gauss.dS_e*param.mesh.d_t*sum(... 234 gauss.alpha_out1(t).*(gauss.Tout1(t) - Te1).*right +... 235 gauss.alpha_out2(t).*(gauss.Tout2(t) - Te2).*right ) ;
236
237 Qin = -gauss.dS_w*param.mesh.d_t*sum(... 238 gauss.alpha_in1(t).*(gauss.Tin1(t) - Tw1).*left +... 239 gauss.alpha_in2(t).*(gauss.Tin2(t) - Tw2).*left ) ;
240
241 %old version 242 %{
243 for e = 1: nel 244 KB1e = zeros (4,4); 245 KB3e = zeros (4,4); 246 RB1e = zeros (4,1); 247 RB3e = zeros (4,1);
```
 %Ce <sup>=</sup> rho( e ) \* gauss . detJ \*(c1( e ).\* gauss . N1\_q +... <sup>249</sup> % c2( e ).\* gauss . N2\_q +c3(e ).\* gauss . N3\_q +... % c4( e ).\* gauss . N4\_q ) ; KCe <sup>=</sup> gauss . detJ \*(k1( e ).\* gauss . B1\_q <sup>+</sup> ... <sup>252</sup> k2( e ).\* gauss . B2\_q <sup>+</sup> k3( e ).\* gauss . B3\_q +... k4( e ).\* gauss . B4\_q ) ; % calculate enthaply if param . mesh . elem\_mat (e ) == <sup>1</sup> <sup>257</sup> H { e } = calc\_h1 (e ) + calc\_h2 (e ) + calc\_h3 (e ) + calc\_h4 (e ) ; H\_l{ e } = hl\_1 ( e ) + hl\_2 ( e ) + hl\_3 ( e ) + hl\_4 ( e ); else <sup>260</sup> H { e } = calc\_h ( e)\*( T1( e ) +T2( e ) +T3( e ) +T4( e ) ) ; end <sup>262</sup> if param . mesh . flag (e ,1) == 1 % right side of PCM <sup>264</sup> nodes <sup>=</sup> param . mesh .IEN (:,e ) ; coord = [ param . mesh . x( nodes ) ; param . mesh . y ( nodes ) ] ; <sup>266</sup> y1 <sup>=</sup> coord (1,2) +( coord (4,2) -coord (1,2)) \*(-0.57735027/2) ; <sup>267</sup> y2 <sup>=</sup> coord (1,2) +( coord (4,2) -coord (1,2)) \*(0.57735027/2) ; <sup>268</sup> KB1e <sup>=</sup> param . icbc . alpha\_out \* gauss . dS\_e \*( ... <sup>270</sup> gauss . Ne1\_q <sup>+</sup> gauss . Ne2\_q ) ; RB1e <sup>=</sup> param . icbc . alpha\_out \* gauss . dS\_e \*... <sup>273</sup> ( param . icbc . Tout (t ,y1) \* gauss .Ne1 +... <sup>274</sup> param . icbc . Tout (t ,y2) \* gauss .Ne2 ) ; % calculate heat flow <sup>277</sup> Qout { e } = gauss . dS\_e \* param . icbc . alpha\_out \* param . mesh .d\_t \*(... param . icbc . Tout (t ,y1) - Te1 ( e ) + ... <sup>279</sup> param . icbc . Tout (t ,y2) - Te2 ( e )) ; elseif param . mesh . flag (e ,1) == 3 % left side of PCM <sup>282</sup> nodes <sup>=</sup> param . mesh .IEN (:,e ) ; coord = [ param . mesh . x( nodes ) ; param . mesh . y ( nodes ) ] ; <sup>284</sup> y1 <sup>=</sup> coord (1,2) +( coord (4,2) -coord (1,2)) \*(-0.57735027/2) ; <sup>285</sup> y2 <sup>=</sup> coord (1,2) +( coord (4,2) -coord (1,2)) \*(0.57735027/2) ; <sup>286</sup> KB3e = -gauss . dS\_w \*( param . icbc . alpha\_in (t ,y1) \* gauss . Nw1\_q ... <sup>288</sup> <sup>+</sup> param . icbc . alpha\_in (t ,y2) \* gauss . Nw2\_q ) ; RB3e = -gauss . dS\_w \*( param . icbc . alpha\_in (t ,y1) \*... <sup>291</sup> param . icbc .Tin(t ,y1) \* gauss .Nw1 +... <sup>292</sup> param . icbc . alpha\_in (t ,y2) \* param . icbc . Tin (t ,y2) ... <sup>293</sup> \* gauss .Nw2 ) ; % calculate heat flow <sup>296</sup> Qin{ e } = -gauss . dS\_w \* param . mesh . d\_t \*(... param . icbc . alpha\_in (t ,y1)\*( param . icbc . Tin (t ,y1) -Tw1 ( e ) )... <sup>298</sup> <sup>+</sup> param . icbc . alpha\_in (t ,y2)\*( param . icbc .Tin(t ,y2) -Tw2( e ) ) ) ; end <sup>300</sup> % convection term <sup>302</sup> if p\_el\_FEM (e ) == <sup>1</sup> <sup>303</sup> KCe <sup>=</sup> KCe <sup>+</sup> u1( e ) \* gauss . NBc1u +v1( e )\* gauss . NBc1v +... <sup>305</sup> u2( e ) \* gauss . NBc2u +v2( e ) \* gauss . NBc2v +... <sup>306</sup> u3( e ) \* gauss . NBc3u +v3( e ) \* gauss . NBc3v +... <sup>307</sup> u4( e ) \* gauss . NBc4u +v4( e ) \* gauss . NBc4v +... <sup>308</sup> u5( e ) \* gauss . NBc5u +v5( e ) \* gauss . NBc5v +... <sup>309</sup> u6( e ) \* gauss . NBc6u +v6( e ) \* gauss . NBc6v +... <sup>310</sup> u7( e ) \* gauss . NBc7u +v7( e ) \* gauss . NBc7v +... <sup>311</sup> u8( e ) \* gauss . NBc8u +v8( e ) \* gauss . NBc8v +... <sup>312</sup> u9( e ) \* gauss . NBc9u +v9( e ) \* gauss . NBc9v ; <sup>313</sup> end <sup>314</sup> K\_e{ e } = reshape (KCe <sup>+</sup> KB1e <sup>+</sup> KB3e ,4\*4,1) ; <sup>315</sup> R\_e{ e } = RB1e <sup>+</sup> RB3e ; <sup>316</sup> C\_e{ e } = reshape (Ce ,4\*4,1) ; <sup>317</sup> end <sup>318</sup> %% assemble to global matrices and vectors <sup>320</sup> % C = sparse ( sp\_mat (:,1) ,sp\_mat (:,2) ,cell2mat (C\_e) ,non ,non) ; % K = sparse ( sp\_mat (:,1) ,sp\_mat (:,2) ,cell2mat (K\_e) ,non ,non) ; % R = sparse ( sp\_vec (:,1) ,sp\_vec (:,2) ,cell2mat (R\_e) ,non ,1); <sup>324</sup> %} end

Code Listing A.8: Code of function c var

```
1 %% c_var 2 function [c]=c_var(T,mat)
  3 % Computes heat capacity c from the value of temperature T 4
  5 c = mat.c_s*( T <= (mat.Tm - mat.dTm) )+... 6 mat.c_l*( T >= (mat.Tm + mat.dTm) )+... 7 (mat.lh + mat.c_s*(mat.Tm + mat.dTm - T) + ... 8 mat.c_l*(T-mat.Tm + mat.dTm))/(2*mat.dTm).*... 9 ( T > (mat.Tm - mat.dTm) & T < (mat.Tm + mat.dTm) );
10 end
```
Code Listing A.9: Code of function k var

```
1 %% k_var
 2 function [k]=k_var(T,mat)
 3 % Computes heat conductivity k from the value of temperature T 4
 5 k = mat.k_s*( T <= (mat.Tm - mat.dTm) )+... 6 mat.k_l*( T >= (mat.Tm + mat.dTm) )+... 7 (mat.k_s*(mat.Tm + mat.dTm - T) + ... 8 mat.k_l*(T-mat.Tm + mat.dTm))/(2*mat.dTm).*... 9 ( T > (mat.Tm - mat.dTm) & T < (mat.Tm + mat.dTm) );
10 end
```
Code Listing A.10: Code of function solve euler back

```
1 %% solve euler back 2 function [Tnew] = solve_euler_back(Told ,C,K,R,param)
 3 % Solves the system of FEM equations via the backwards euler method 4
 5 Tnew = Told +param.mesh.d_t*... 6 ( (C + param.mesh.d_t*K)\(R-K*Told) );
 7 % long version : 8 %{
9 % build equation system A * X = b 10 M_eff = C + param . mesh .d_t* K ;
11 b = R - K * Told (:); 12
13 % solve equation system A * X = b 14 X = M_eff \ b ;
15 Tnew = Told + param . mesh .d_t * X;
16 %}
```
Code Listing A.11: Code of function waitbar2D simple

```
1 %% waitbar 2 D simple
 2 function [f,timearray ,up] = waitbar2D_simple(f,param ,timearray ,toc ,t,up)
 3 % updates waitbar in a simple way
 4
 5 wait = t/param.mesh.not;
 6 timearray(t) = toc -sum(timearray);
 7 % moving average of computation time 8 average = movmean(timearray ,[20 0], Endpoints ,timearray( end));
 9 waittime = (param.mesh.not -t)*average(end );
10
11 if (toc -up) > param.display.update
12 % wait at least a specified time between updates
13 up = toc;
14
15 if waittime < 60 16 waitbar(wait ,f,sprintf(... 17 Fortschritt : %2.1 f %%, Restdauer : ~ %2.1 f s ,100*wait ,waittime));
18 else 19 waittime_s = mod(waittime ,60);
20 waittime_m = floor(waittime/60);
21 waitbar(wait ,f,sprintf(... 22 Fortschritt : %2.1 f %%, Restdauer : ~ %3.0 f min %2.0 f s ,... 23 100*wait ,waittime_m ,waittime_s));
24 end 25 end 26 end
```
Code Listing A.12: Code of function solve navierstokes

```
1 %% solve navier stokes 2 function [U_new ,V_new ,grad ,U,V] = solve_navierstokes(U,V,T_ij ,p_el ,param)
 3 % Solves the navier stokes equations via finite difference method 4
 5 %% preprocessing
 6 nx = param.mesh.nx;
 7 ny = param.mesh.ny;
 8 f_u = reshape ((p_el (1:nx -1,:) & p_el (2:nx ,:)) ,[],1);
 9 f_v = reshape ((p_el (:,1:ny -1) & p_el (:,2:ny)) ,[],1);
10 % temperatures for boussinesque approximation
11 T_v = reshape(avg(T_ij(:,2: end -1)) ,[],1);
12 T_u = reshape(avg(T_ij(2: end -1,:) ) ,[],1);
13
14 f2_v = 0+(T_v >=param.pcm.Tm+param.pcm.dTm).*1+... 15 (T_v >param.pcm.Tm -param.pcm.dTm & T_v <param.pcm.Tm+param.pcm.dTm).*... 16 (T_v -param.pcm.Tm+param.pcm.dTm)/(2* param.pcm.dTm);
17 f2_u = 0+(T_u >=param.pcm.Tm+param.pcm.dTm).*1+... 18 (T_u >param.pcm.Tm -param.pcm.dTm & T_u <param.pcm.Tm+param.pcm.dTm).*... 19 (T_u -param.pcm.Tm+param.pcm.dTm)/(2* param.pcm.dTm);
20
21 %% nonlinear terms / explicit time step
22 gamma = min(1.2*param.mesh.d_t*max( ... 23 max(max(abs(U)))/param.geom.dx,max(max(abs(V)))/param.geom.dy ) ,1);
24 % add boundary conditions 25 Ue = [zeros(1,ny+2); [-U(:,1) U -U(:, end)] ;zeros(1,ny+2)];
26 Ve = [0 -V(1,:) 0; [zeros(nx ,1) V zeros(nx ,1)];0 -V(end ,:) 0];
27
28 % upwinding
29 Ua = avg(Ue ) ; % vertical average of U 30 Ud = diff(Ue ) /2; % vertical finite difference of U 31 Va = avg(Ve); % horizontal average of V 32 Vd = diff(Ve)/2; % horizontal finite difference of V 33 UVx = diff(Ua.*Va -gamma*abs(Ua).*Vd)/param.geom.dx; % d (UV) /dx 34 UVy = diff((Ua.*Va -gamma*Ud.*abs(Va)) ) / param.geom.dy; % d (UV) /dy
35 Ua = avg(Ue(:,2: end -1)); % horizontal average of U 36 Ud = diff(Ue(:,2: end -1))/2; % horizontal finite difference of U 37 Va = avg(Ve(2: end -1,:) ) ; % vertical average of V 38 Vd = diff(Ve(2: end -1,:) ) /2; % vertical finite difference of V 39 U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/param.geom.dx; % d (UU) /dx 40 V2y = diff((Va.^2-gamma*abs(Va).*Vd) ) / param.geom.dy; % d (VV) /dx 41 % explicit time step
42 U = U-param.mesh.d_t*(UVy(2: end -1,:)+U2x);
43 V = V-param.mesh.d_t*(UVx(:,2: end -1)+V2y);
44
45 %% viscosity terms / implicit time step
46 % implicit step for U 47 Ex=sparse (2:nx -1,1:nx -2,1,nx -1,nx -1);
48 Ax=Ex+Ex -2* speye(nx -1);
49 Ey=sparse (2:ny ,1:ny -1,1,ny,ny);
50 Ay=Ey+Ey -2* speye(ny);
51 Ay(1,1)=-3; Ay(ny,ny)=-3; % Neumann B .Cs 52 A=kron(Ay/param.geom.dy^2,speye(nx -1))+kron(speye(ny),Ax/param.geom.dx^2);
53
54 Uo=reshape(U,[],1);
55 mue = param.pcm.mue*(1+(1-f2_u).^4*10^3);
56 Fu = (1-f_u).^2./( f_u .^3+0.001) *10^8;
57 Un = (speye(length(Uo)).*(1+Fu) - param.mesh.d_t/param.pcm.rho*mue.*A) \ ... 58 (Uo - sin(param.pcm.phi)*param.mesh.d_t*param.pcm.grav*... 59 param.pcm.beta*(param.pcm.T_0 - T_u));
60 U=reshape(Un,nx -1,ny);
61
62 % implicit step for V 63 Ex=sparse (2:nx ,1:nx -1,1,nx,nx);
64 Ax=Ex+Ex -2* speye(nx);
65 Ax(1,1)=-3; Ax(nx,nx)=-3; % Neumann B .Cs 66 Ey=sparse (2:ny -1,1:ny -2,1,ny -1,ny -1);
67 Ay=Ey+Ey -2* speye(ny -1);
68 A=kron(Ay/param.geom.dy^2,speye(nx))+kron(speye(ny -1),Ax/param.geom.dx^2);
69
70 Vo=reshape(V,[],1);
71 mue = param.pcm.mue*(1+(1-f2_v).^4*10^3);
72 Fv = (1-f_v).^2./( f_v .^3+0.001) *10^8;
73 Vn = (speye(length(Vo)).*(1+Fv) - param.mesh.d_t/param.pcm.rho*mue.*A) \ ... 74 (Vo - cos(param.pcm.phi)*param.mesh.d_t*param.pcm.grav*... 75 param.pcm.beta*(param.pcm.T_0 - T_v));
76 V=reshape(Vn ,nx ,ny -1);
77
```

```
78 %% pressure correction 79 grad = reshape(diff([zeros(1,ny);U;zeros(1,ny)])/param.geom.dx... 80 +diff([zeros(nx ,1) V zeros(nx ,1) ] ) / param.geom.dy ,[],1);
81
82 Lp_p = Laplace_p(p_el ,param);
83 p = Lp_p\(param.pcm.rho/param.mesh.d_t*grad);
84
85 P = reshape(p,nx,ny);
86 U_new = U-param.mesh.d_t/param.pcm.rho*diff(P)/param.geom.dx;
87 V_new = V-param.mesh.d_t/param.pcm.rho*diff(P ) / param.geom.dy;
88
89 end
```
Code Listing A.13: Code of function Laplace p

```
1 %% Laplace p
2 function [Lp_p] = Laplace_p(p_el ,param)
3 % Computes the finite difference 2 D Laplace stencil for the pressure p .
4 % Includes homogeneous Neumann boundary conditions on the internal phase
5 % boundaries . 6
7 dxq = param.geom.dx^2;
8 dyq = param.geom.dy^2;
9 nx = param.mesh.nx;
10 ny = param.mesh.ny;
11
12 ps = reshape ([false(nx ,1) p_el (:,1: end -1)],[],1);
13 pn = reshape([p_el(:,2: end) false(nx ,1)],[],1);
14 pw = reshape([false(1,ny); p_el(1: end -1,:)],[],1);
15 pe = reshape([p_el(2: end ,:); false(1,ny)],[],1);
16 p_el = reshape(p_el ,[],1);
17
18 i = param.mesh.stag(:,1);
19 j = param.mesh.stag(:,2);
20
21 Lp_val = [i+(j-1)*nx, i+(j-1)*nx,... 22 0.1-1/dxq*(p_el==pw) -1/dxq*(p_el==pe) -1/dyq*(p_el==ps) -1/dyq*(p_el
                   ==pn);
23 i+(j-1)*nx, i-1+(j-1)*nx, 0+1/ dxq*(p_el==pw);
24 i+(j-1)*nx, i+1+(j-1)*nx, 0+1/ dxq*(p_el==pe);
25 i+(j-1)*nx, i+(j-2)*nx, 0+1/ dyq*(p_el==ps);
26 i+(j-1)*nx , i+(j)*nx , 0+1/ dyq*(p_el==pn);
27 ];
28 % small Value 0.1, so that Lp_p is not singular . 29 Lp_val = Lp_val(find(( Lp_val (:,2) >0)&( Lp_val (:,2) <=param.mesh.nel)) ,:);
30
31 Lp_p = sparse(Lp_val(:,1),Lp_val(:,2),Lp_val(:,3),param.mesh.nel ,param.mesh.
         nel);
32
33 end
```
# Appendix B Additional plots

To illustrate the simulation results of the parameter study discussed in Section 6.2, additional plots which show the temperature and velocity distribution in the simulated cavities at six time points during the respective simulations are presented here.

For the sake of improving comparability of the given plots, the illustrations begin on the following double page, so that each double page contains the temperature and velocity distributions of a particular cavity orientation case.

Figure B.1: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 0

Figure B.1: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 0

Figure B.2: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 45

Figure B.2: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 45

Figure B.3: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 90

Figure B.3: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 90

Figure B.4: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 135

Figure B.4: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 135

Figure B.5: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 180

Figure B.5: Temperature and velocity distribution at di↵erent simulation times for cavity orientation 180

To increase the effciency of energyintensive industrial processes, thermal energy storages can offer new possibilities. A novel approach is investigated in the project HyStEPs. In this concept, containers flled with PCM are placed at the shell surface of a Ruths steam storage, to increase storage effciency.

In this work, a two-dimensional model using the fnite element method is developed to simulate the PCM of the hybrid storage as designed in the HyStEPs project. The apparent heat capacity method is applied in a MATLAB implementation, considering heat transfer by both conduction and natural convection. This successfully validated code can handle any desired layout of materials arranged on a rectangular domain.

Kasper

**HYBRID STORAGE USING THE FINITE ELEMENT METHOD**

**MODELING OF THE PHASE CHANGE MATERIAL OF A** 

Furthermore, a parameter study of different dimensions and orientations of the PCM cavity was conducted. The impact of natural convection was found to lead to signifcantly varying behaviour of the studied cavities with different orientation during the charging process, while it was found to be negligible during the discharging process.