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
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.
Hy
StE
Ps
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
[email protected]
www.tuwien.at/academicpress
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License
(CC BY-SA 4.0). https://creativecommons.org/licenses/by-sa/4.0/
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
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.
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 Efficient Processes [HyStEPs]. The aim
of HyStEPs is to develop an all-in-one storage concept to increase the efficiency 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 con-
siderable 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.
v
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. René Hofmann
January, 2020
vi
Acknowledgements
I am very grateful to my thesis advisor, Prof. René 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.
vii
viii
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 efficiency 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 efficiency 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 rectan-
gular 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 beha-
viour 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 vary-
ing 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.
ix
x
Table of Contents
Nomenclature xiii
List of Figures xvii
List of Tables xxi
List of code Listings xxiii
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Aim of this work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Methodological approach . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5 Structure of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Theoretical framework 7
2.1 Basics of thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Discretization methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Numerical modelling 21
3.1 Model equations and assumptions . . . . . . . . . . . . . . . . . . . . . 22
3.2 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 Convection modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
xi
Table of Contents
4 MATLAB implementation 39
4.1 Structure of the program . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2 Description of the key functions . . . . . . . . . . . . . . . . . . . . . . 42
4.3 Usage information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5 Validation 51
5.1 Validation of heat conduction model . . . . . . . . . . . . . . . . . . . 52
5.2 Validation of convection model . . . . . . . . . . . . . . . . . . . . . . 57
5.3 Cross-validation with ANSYS Fluent . . . . . . . . . . . . . . . . . . . 64
5.4 Test of a fluctuating charging/discharging scenario . . . . . . . . . . . 67
6 Parameter studies 71
6.1 Study of di↵erent cavity designs . . . . . . . . . . . . . . . . . . . . . . 71
6.2 Comparison of di↵erent cavity orientations . . . . . . . . . . . . . . . . 78
6.3 Discussion of conduction vs. convection simulations . . . . . . . . . . . 82
7 Conclusion 83
7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.2 Scope and limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.3 Suggested future objectives . . . . . . . . . . . . . . . . . . . . . . . . 87
Bibliography 91
A MATLAB code listings 97
B Additional plots 113
xii
Nomenclature
Variables and Parameters
⇥ ⇤
↵ heat transfer coefficient W/(m2 K)
⇥ ⇤
↵L thermal di↵usivity m2 /s
linear thermal expansion coefficient [1/K]
lm specific latent heat [J/kg]
t simulation time step size [s]
x, y mesh resolution in x and y direction [m]
✏ mushy region temperature range [K]
St dimensionless parameter occuring in the solution of the Stefan
problem
⇥ ⇤
f force density vector N/m3
⇥ ⇤
g gravitational acceleration vector m/s2
u velocity vector [m/s]
⇥ ⇤
µ dynamic viscosity (N s)/m2
⇥ ⇤
⇢ density kg/m3
⇥ ⇤
⇢0 reference density kg/m3
⇠, ⌘ transformed space coordinates [m]
am fraction of total mass
ar aspect ratio
c specific heat capacity [J/(kgK)]
capp apparent heat capacity [J/(kgK)]
cp specific heat capacity at constant pressure [J/(kgK)]
H enthalpy [J]
k thermal conductivity [W/(mK)]
m mass [kg]
n number of simulation/measurement values
⇥ ⇤
p pressure N/m2
Q heat storage capacity [J]
xiii
Nomenclature
⇥ ⇤
q heat flux W/m2
s melting front location [m]
StL Stefan number
T temperature [K]
t time [s]
u velocity component in x direction [m/s]
v velocity component in y direction [m/s]
x, y space coordinates [m]
Xaluminium cavity aluminium containment width [m]
XP CM cavity PCM dimension in x direction [m]
Xsteel steel containment width [m]
Yaluminium cavity aluminium fin width [m]
YP CM cavity PCM dimension in y direction [m]
B derivative of shape functions N
N shape functions
Subscripts
0 initial
boundary value at the boundary
el number of elements
in input value
liquid material parameter of liquid phase
m melting point
max maximum value
n number of nodes
out output value
ref reference value
solid material parameter of solid phase
Superscripts
⇤ ⇤⇤
, intermediate solutions of the fractional step approach
h horizontal averaged/di↵erenced quantity
v vertical averaged/di↵erenced quantity
xiv
Nomenclature
Abbreviations
CFD computational fluid dynamics
CFL Courant-Friedrichs-Lewy (condition)
FDM finite di↵erence method
FEM finite element method
FVM finite volume method
HyStEPs hybrid storage for efficient processes
PCM phase change material
RMSE root mean square error
RSS Ruths steam storage
Symbols
finite di↵erence/Laplace operator = r2
r nabla operator: r = (@/@x, @/@y)
xv
xvi
List of Figures
1.1 Schematic sketch of the hybrid storage design consisting of the Ruths
steam storage and the surrounding PCM modules . . . . . . . . . . . . 3
2.1 Enthalpy H versus temperature T with (a) mushy phase change and
(b) isothermal phase change . . . . . . . . . . . . . . . . . . . . . . . . 9
3.1 Four-noded bilinear element in natural coordinate system ⇠, ⌘ with

node coordinates in brackets . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Two-dimensional domain discretized into a staggered grid . . . . . . . 35
4.1 Flowchart diagram of FEM model PCMsolver . . . . . . . . . . . . . . 41
4.2 Image of a 2x2 mesh, with numbered nodes and elements . . . . . . . 43
5.1 Time evolution of the melting front of simulations with di↵erent val-
ues of the mushy region size ✏ compared to the melting front of the
analytical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.2 Time evolution of liquid fraction for di↵erent aspect ratio values ar =
x
relative to the base case x = y = 0.5 mm compared to experi-
y
mental data with variable x . . . . . . . . . . . . . . . . . . . . . . . 61
5.3 Time evolution of liquid fraction for di↵erent aspect ratio values ar =
x
relative to the base case x = y = 0.5 mm compared to experi-
y
mental data with variable y . . . . . . . . . . . . . . . . . . . . . . . 62
5.4 Melting fronts of two simulations with di↵erent time step size compared
to experimental data, given at five time points . . . . . . . . . . . . . 63
5.5 Melting fronts of two simulations with di↵erent mesh size compared to
experimental data, given at five time points . . . . . . . . . . . . . . . 64
5.6 Phase structure in PCM cavity at di↵erent times . . . . . . . . . . . . 66
xvii
List of Figures
5.7 Temperature and velocity evolution over time in fluctuating gallium

test case charging/discharging scenario . . . . . . . . . . . . . . . . . . 68
5.8 Temperature distribution in the gallium cavity at the simulation time

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

minutes during charging/discharging scenario . . . . . . . . . . . . . . 69
6.1 Sketch of a typical PCM cavity fin section attached to the steel con-
tainment of the steam storage . . . . . . . . . . . . . . . . . . . . . . . 72
6.2 Enthalpy per surface area of the simulated cases with fin width
Yaluminium = 5 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.3 Enthalpy per surface area of the simulated cases with fin width
Yaluminium = 2.5 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.4 Enthalpy per surface area of the simulated cases with di↵erent fin width
Yaluminium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.5 Schematic illustration of the five simulated cavity orientations in rela-

tion to the gravity vector . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.6 Change of enthalpy during the charging process . . . . . . . . . . . . . 81
6.7 Evolution of the melted volume fraction during charging (a) and dis-
charging (b) of the cavity for di↵erent orientations . . . . . . . . . . . 82
B.1 Temperature and velocity distribution at di↵erent simulation times for

cavity orientation 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

cavity orientation 45 . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

xviii
List of Figures


xix
xx
List of Tables
3.1 Thermophysical material properties . . . . . . . . . . . . . . . . . . . . 23
3.2 Heat transfer coefficients ↵ occurring between the inside of the RSS
and its outer wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3 Gauss points and weight factors of the used quadrature rules . . . . . 31
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) . . . . . . . . . . . . . . . . . . . . . 53
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) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.3 RMSE of melting front position for di↵erent values of mushy region
size ✏ compared to Stefan solution (analytical) . . . . . . . . . . . . . 54
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) . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.5 RMSE of liquid fraction for di↵erent values of the time step size com-
pared to test case (experimental) and t = 2 ms simulation (convergence) 59
5.6 RMSE of liquid fraction for di↵erent values of the mushy region size
compared to test case (experimental) . . . . . . . . . . . . . . . . . . . 60
6.1 Simulated parameter study cases and associated tags with fin width
Yaluminium = 5 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.2 Simulated parameter study cases and associated tags with fin width
Yaluminium = 2.5 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
xxi
xxii
List of code Listings
4.1 Listing of a typical input datafile input datafile.m . . . . . . . . . . . . 46
A.1 Code of main function PCMsolver . . . . . . . . . . . . . . . . . . . . . 97
A.2 Code of function mesh2d . . . . . . . . . . . . . . . . . . . . . . . . . . 99
A.3 Code of function shape mat . . . . . . . . . . . . . . . . . . . . . . . . 101
A.4 Code of shape matrix function NmatHeat2D . . . . . . . . . . . . . . . 102
A.5 Code of shape matrix function BmatHeat2D . . . . . . . . . . . . . . . 102
A.6 Code of FEM preprocessing function FEM preprocessor . . . . . . . . . 102
A.7 Code of FEM matrix assembly function heat2Delem new . . . . . . . . 104
A.8 Code of function c var . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
A.9 Code of function k var . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
A.10 Code of function solve euler back . . . . . . . . . . . . . . . . . . . . . 110
A.11 Code of function waitbar2D simple . . . . . . . . . . . . . . . . . . . . 110
A.12 Code of function solve navierstokes . . . . . . . . . . . . . . . . . . . . 111
A.13 Code of function Laplace p . . . . . . . . . . . . . . . . . . . . . . . . . 112
xxiii
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 efficiency 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 efficient 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 System-
technik GmbH.
The basic approach of HyStEPs is to encase Ruths steam storages in innovative lat-
ent heat accumulators. Furthermore, a state of charge measurement for the hybrid
storage tank should ensure optimal control and operation of both storage types dur-
ing charging and discharging. This innovative solution is based on an AIT patent
application.
1
Chapter 1. Introduction
1.2 Problem statement
Recently, the applications for latent heat thermal energy storages have increased sub-
stantially. 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 difficult 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.
2
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
requirements. In addition, the vessels manufactured must fulfil the requirements
Aim offor
thisapprovals,
work tests,
examinations, and inspections.
PCM container
Steam
Liquid
Figure 6. Schematic sketch of the hybrid storage design consisting of the RSS and the surrounding
Figure
LHTES 1.1: [43].
modules Schematic sketch of the hybrid storage design consisting of the Ruths
steam storage and the surrounding PCM modules [2]
1.3 Aim of this work
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
3
Figure 7. Schematic sketch of the PCM container mounting: (a) front view. (b) top view. (c) sectional
view A-A. According to [43].
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 char-
ging/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, modifica-
tions 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, approx-
imation 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 temper-
ature 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 efficient matrix manipulation
operations. This is especially suitable for the finite element method.
4
Structure of this thesis
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 para-
maters 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
5
the comparative data was obtained and the developed model was therefore success-
fully 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 charac-
teristics 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.
6
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
Z Tf
Q= mcp dT (2.1)
Ti
depends on the specific heat capacity cp at constant pressure p of the material, the
temperature di↵erence (Tf 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
Z Tm Z Tf
Q= mcp,solid dT + mam lm + am mcp,liquid dT , (2.2)
Ti Tm
with the specific latent heat lm and the liquid fraction of the total mass am [6]. The
specific heat capacity at constant pressure cp 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
7
Chapter 2. Theoretical framework
latent heat and small volume changes and o↵er the advantages of less stringent
container requirements and greater design flexibility [7].
Solid-gas transitions, as well as direct liquid-gas transitions, have very high

latent heat of phase transition, but their large volumetric expansion during
these transitions is associated with containment problems and often rule out
their potential utility in thermal energy storage systems [6].
Solid-liquid transformations have comparatively smaller latent heat than liquid-

gas transformations. However, these transformations involve only a small change
(on the order of 10% or less) in volume and have therefore proven to be eco-
nomically attractive for use in thermal energy storage systems [6].
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 Tm [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
8
Basics of thermodynamics
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
q= krT , (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
q = ⇢cT u , (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
@(⇢cT )
= r(krT ) r(⇢cT u) . (2.5)
@t
9
According to the literature, this equation is also called the convection-di↵usion equa-
tion 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.
@(⇢)
= r(⇢u) (2.6)
@t
@(⇢)
For an incompressible fluid with @t = 0, this continuity equation takes the form
r(⇢u) = 0 . (2.7)
The momentum balance, also known as the Navier-Stokes-equation, takes the fol-
lowing form for an incompressible Newtonian fluid with constant material properties
µ = const., ⇢ = ⇢0 = const. :
✓ ◆
@u
⇢0 + (u r)u r (µru) + rp = f [13] , (2.8)
@t
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 ( ⌦, t) = Tboundary (t) . (2.9)
Neumann boundary conditions prescribe a heat flux qboundary (t) at the bound-
ary. By using Fourier’s law (2.3) these can be written as
@T
k | ⌦,t = qboundary (t) , (2.10)
@n
10
Basics of thermodynamics
where n is the outward unit normal to the boundary surface [9]. For perfectly
insulated boundaries (adiabatic wall) this means
@T
k | ⌦,t = 0. (2.11)
@n
Robin boundary conditions prescribe a heat exchange with the surrounding

boundary at the ambient temperature Tboundary (t) with a specified thermal con-
ductivity ↵boundary :
krT | ⌦,t = ↵boundary (T | ⌦,t Tboundary (t)) . (2.12)
This type of boundary condition is often also called Fourier condition or con-
vective 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
u| ⌦,t = 0, (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.
11
The simplest form of a mathematical model describing phase transitions is the clas-
sic, one-dimensional Stefan problem. This problem describes the liquefaction of a
solid medium, initially at the solidification temperature T (x, t = 0) = Tm with the
boundary condition T (x = 0, t) = T0 > 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))
= r(krT (x, t)) heat equation, 0 < x < s(t), t > 0
@t
T (x = 0, t) = T0 > Tm Boundary condition, t > 0
T (x, t = 0) = Tm Initial condition
The free boundary x = s(t)

@s @T
lm ⇢ = k Stefan condition
@t @x
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
c(T0 Tm ) k
number StL = lm and the thermal di↵usivity ↵L = c⇢
T0 x
T (x, t) = T0 erf ( p ) (2.14)
erf ( St ) 2 ↵L t
p
s(t) = 2 St ↵L t (2.15)
2 StL
St e
St erf ( St ) = p . (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.
12
Numerical methods
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:
ZT
H(T ) = H(T ) + ⇢f (T )a = ⇢c(T )dT + ⇢f (T )a (2.17)
Tref
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
8
>
> 0 , solid: T < Tm
>
<
a= ] 0, 1 [ , mushy: T = Tm (2.18)
>
>
>
: 1 , liquid: T > T m
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)
✓ ◆
@ H @ @ H @f
= ↵L ⇢a , (2.19)
@t @x @x @t
k
with the thermal di↵usivity ↵L = c⇢ .
The enthalpy method leads to a simple phase change problem, since interface condi-
tions must not be taken into account [17].
13
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:
8
>
> c , solid: Tm " > T
>
< solid
lm +csolid (Tm +" T )+cliquid ( Tm +"+T )
capp = , mushy: Tm "  T  Tm + "
>
> 2"
>
: cliquid , liquid: T > Tm + " .
(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 meth-
ods 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
14
Discretization methods
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 con-
tinuous 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 ] {T } . (2.21)
The general approach of the finite element method is to multiply the partial di↵er-
ential equation (PDE) with given boundary conditions (the strong form) by a test
function, or weighting function, and integrate over the whole simulation domain. Ap-
plying Green’s first integration theorem (partial integration) leads to a variational
formulation, also called the weak form [20]. To discretize this still infinite dimen-
sional 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 estab-
lished in computational fluid dynamics (CFD) applications, including the well-known
15
software package ANSYS Fluent, Nedjar [9] states that especially finite element meth-
ods (FEM) are able to handle complex coupled thermomechanical problems with vari-
ous 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]
@u u(x + x) u(x) u(x) u(x x)

(x) = lim = lim =
@x x!0 x x!0 x
u(x + x) u(x x)
lim (2.22)
x!0 2 x
for an infinitesimal di↵erence x between the grid points. This leads to three types
of commonly used approximations:

@u ui+1 ui
⇡ forward di↵erence
@x x
 i
@u ui ui 1
⇡ backward di↵erence (2.23)
@x x
 i
@u ui+1 ui 1
⇡ central di↵erence .
@x i 2 x
The derivatives, which typically arise in two-dimensional convection-di↵usion prob-

lems 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:

@u ui+1,j ui 1,j
⇡ (2.24)
@x i,j 2 x

@u ui,j+1 ui,j 1
⇡ (2.25)
@y i,j 2 y

@2u ui+1,j 2ui,j + ui 1,j
⇡ (2.26)
@x2 i,j 4 x2

@2u ui+1,j+1 ui+1,j 1 ui 1,j+1 + ui 1,j 1
⇡ . (2.27)
@x@y i,j 4 x y
16
The derivative of a nonlinear term u(x, y) v(x, y) can be calculated via the product
rule thus
  
@uv @v @u
= ui,j + vi,j ⇡
@y i,j @y i,j @y i,j
ui,j vi,j
(vi,j+1 vi,j 1) + (ui+1,j ui 1,j ) . (2.28)
2 y 2 x
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 compu-
tational 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 tak-
ing the integral form of conservation equations and splitting the integrals into small
intervals. The finite volume method is strictly conservative, because, since two neigh-
bouring 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
@T T n+1 T n
= . (2.29)
@t t
Here, T n = T (x, tn ) is the variable’s value at t = tn and t is the time increment
according to tn+1 = tn + t. The variable T (x, tn ) is usually assumed to be known
17
and used as an initial condition to advance the solution to the next time level. The
solution T can then be written as
T = ✓T n+1 + (1 ✓)T n , tn  t  tn+1 . (2.30)
The relaxation parameter ✓ 2 [0, 1] is most commonly one of the following:
✓=1 implicit backwards Euler scheme

1
✓= 2 second-order, centered implicit method (Crank-Nicolson-Galerkin) [19]
✓=0 explicit forward Euler scheme
(2.31)
1
It is known that for ✓ < 2, 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

n o
[C] Ṫ + [K] {T } = [R] (3.23)
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 :
n o
[C] Ṫ + [K] {T }n+1 = [R] . (2.32)
n+1
The time derivative of the temperature vector {T } is approximated through the back-
wards di↵erential quotient
⇢ n o
@T {T }n+1 {T }n
= Ṫ = , (2.33)
@t n+1 n+1 t
which gives us the value

n o
{T }n+1 = {T }n + t Ṫ . (2.34)
n+1
Combining equation (2.32) with (2.34) leads, with the expression
[M ] = [C] + t [K] , (2.35)
18
to
n o n o
[M ] Ṫ = ([C] + t [K]) Ṫ = [R] [K] {T }n . (2.36)
n+1 n+1
n o
We can therefore derive the value for Ṫ as
n+1
n o
1
Ṫ = [M ] ([R] [K] {T }n ) , (2.37)
n+1
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
n o
{T }n+1 = {T }n + t Ṫ
n+1
{T }n + t [M ]
1
([R] [K] {T }n ) = (2.38)
1
{T }n + t ([C] + t [K]) ([R] [K] {T }n ) .
19
20
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 condi-
tions 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 simu-
late 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.
21
Chapter 3. Numerical modelling
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 geo-
metry and operating temperatures of this hybrid storage system is a eutectic mixture
of potassium nitrate and sodium nitrate (KNO3 -NaNO3 ). The properties of this mix-
ture, 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 Tm = 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
22
Model equations and assumptions
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].
Table 3.1: Thermophysical material properties

Note. (s) denotes the solid phase, (l) denotes the liquid phase.
⇢ c k Tm lm µ
Material 3 1 1 1 1 2
kg m J(kg K) W (m K) C kJ kg K N sm
KNO3 - 2050(s) 1350(s) 0.457(s) 4 3

220 108 3.5 10 5.8 10
NaNO3 1959(l) 1492(l) 0.435(l)
Steel 7800 540 51 - - - -
Aluminium 2700 910 237 - - - -
Air 1 1000 0.024 - - - -
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 difficulties 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 coefficient 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 coefficients are listed in Table 3.2.
23
In the parameter studies following in Chapter 6, heat transfer coefficients were as-
sumed that are within the range given here. However, the impact of di↵erent heat
transfer coefficients 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 coefficients ↵ occurring between the inside of the RSS and
its outer wall [32]
↵
Operating mode Phase Dominating heat transfer e↵ect 2 1
Wm K
charging vapour condensing steam 5000
charging liquid liquid water 700
discharging vapour dry steam 10
discharging liquid evaporating liquid water 1000
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:
@(⇢cT )
energy equation: = r(krT ) r(⇢cT u) (2.5)
@t
continuity equation: r(⇢u) = 0 (2.7)
✓ ◆
@u
Navier-Stokes equation: ⇢ + (u r)u r (µru) + rp = f (2.8)
@t
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
@(T )
⇢c = r( q) ⇢cr(T u) , (3.1)
@t
24
Model equations and assumptions
with the flux q = krT . Assuming constant density also means freedom of diver-
gence of the velocity field, thus equation (2.7) is reduced to
ru = 0 . (3.2)
To account for phase change in the material, the apparent heat capacity method is
used, as introduced in Section 2.2.2:
8
>
> c , solid: Tm ">T
>
< solid
lm +csolid (Tm +" T )+cliquid ( Tm +"+T )
capp = , mushy: Tm "  T  Tm + "
>
> 2"
>
: c
liquid , liquid: T > Tm + "
(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 particu-
lar form, that of an incompressible Newtonian fluid with constant material properties
⇢, µ = const. A further approximation is to assume constant pressure and only ac-
count for density di↵erences in terms of the force density f . This is known as the
Boussinesq approximation, where the buoyancy force
f= ⇢g = ⇢0 ( (T Tref )) g (3.3)
can be calculated by a linear thermal expansion coefficient , the reference density ⇢0

and a reference temperature Tref [33]. Thus, the Navier-Stokes equation to be solved
is ✓ ◆
@u
⇢0 + (u r)u µr (r u) + rp = f . (3.4)
@t
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:
25
Robin boundary conditions (see (2.12)) have been applied to the east and west
domain boundaries:
|q(x = 0, y)| = |qin (y)| = ↵in (T (x) Tin ) and (3.5)
|q(x = Lx , y)| = |qout (y)| = ↵out (T (x) Tout ) . (3.6)
Adiabatic Neumann boundary conditions have been applied to the south and
north domain boundaries:
q(x, y = 0) = 0 and (3.7)
q(x, y = Ly ) = 0 . (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 Lx and Ly
denote the length of the domain in the corresponding dimension.
Regarding the velocity field u = [u, v]T , no slip boundary conditions are set for the
domain boundaries
u = [u, v]T = 0 , for x = Lx , y = Ly , x = 0, y = 0 , (3.9)
and the velocity is also set to zero if the PCM is solid, as well as for every non-melting
material
u = [u, v]T = 0 , for T < Tm ✏. (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 discret-
ization, which will be explained in Section 3.3.
26
Discretization
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 = nx ny finite elements. In this thesis, four-noded bilinear rectangular
elements are used, which leads to the total sum of nnon = (nx + 1) (ny + 1). The
dependent variables, meaning the temperature T and the velocity u = [u, v]T , can
therefore be approximated by their nodal values as
T = [N ] {T } , (3.11)
u = [N ] {u} and (3.12)
v = [N ] {v} , (3.13)
where [N ] = [N1 , N2 , N3 , ..., 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
q= k [B] {T } (3.14)
when using the gradient of the shape function matrix

2 @N 1 @N 2 @N 3 3
...
6 @x @x @x 7
[B] = 6
4 @N 1
7 .
5 (3.15)
@N 2 @N 3
...
@y @y @y
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
(⇠, ⌘) = N1 1 + N2 2 + N3 3 + N4 4 , (3.16)
27
with the shape functions
1
N1 = (1 ⇠)(1 ⌘) ,
4
1
N2 = (1 + ⇠)(1 ⌘) ,
4 (3.17)
1
N3 = (1 + ⇠)(1 + ⌘) ,
4
1
N4 = (1 ⇠)(1 + ⌘) [19] .
4
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

Z Z Z
T @T T T
[N ] ⇢c dV + [N ] r(q) dV + [N ] ⇢cr(T u) dV = 0 . (3.18)
V @t V V
R R
When applying the divergence theorem V rF dV = S F n dS [34] on the third
integral we obtain the following relation:
Z Z Z
T T T
[N ] ⇢cr(T u) dV + r([N ] )⇢cT u dV = ⇢cr([N ] T u) dV =
V V
ZV (3.19)
T
⇢c([N ] T u)T n dS .
S
28
Discretization
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:
Z Z
T T
[N ] ⇢cr(T u) dV = r([N ] )⇢cT u dV . (3.20)
V V
Similarly, inserting the flux q = krT into the second integral and applying the
divergence theorem leads to
Z Z Z
T T T
[N ] r(q) dV = r([N ] )krT dV + [N ] q n dS , (3.21)
V V S
which is written as
Z Z
T T
[B] k [B] {T } dV + [N ] ↵out ([N ] {T } Tout ) dS+
V S1
Z Z Z (3.22)
T
[N ] ↵in (Tin [N ] {T }) dS + 0 dS + 0 dS
S3 S2 S4
when applying the boundary conditions for the boundary S, where the index S1
denotes the east boundary, S3 the west boundary and indices S2, S4 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
n o
[C] Ṫ + [K] {T } = [R] , (3.23)
with the following system matrices:

Z
T
[C] = [N ] ⇢c [N ] dV , (3.24)
V
Z Z
T T
[K] = [B] k [B] dV ([B] )⇢cu [N ] dV
V V
Z Z
T T
+ ↵out ([N ] [N ] dS ↵in ([N ] [N ] dS , (3.25)
S1 S3
Z Z
T T
[R] = ↵out Tout [N ] dS ↵in Tin [N ] dS . (3.26)
S1 S3
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]T and the thermal conductivity
k = [kx , ky ]T = k [1, 1]T are evaluated. The time discretization of equation (3.23)
29
is done by application of the implicit backwards Euler method as derived in Section

2.3.4 of this thesis:
n o
{T }n+1 = {T }n + t Ṫ
n+1
{T }n + t [M ]
1
([R] [K] {T }n ) = (2.38)
1
{T }n + t ([C] + t [K]) ([R] [K] {T }n ) .
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:
Z Z 1 Z 1
f (x, y) dV = f (⇠, ⌘) |J| d⇠d⌘ (3.27)
V 1 1
The factor |J| is the determinant of the Jacobian matrix of the transformation J [19],
which is defined by
2 @N 3 2 @x @y 3 2 @N 3 2 @N 3
6 @⇠ 7 6 @⇠ @⇠ 7 6 @x 7 6 @x 7
6 7 6 7 6 7=J 6 7 .
4 @N 5 = 4 @x @y 5 4 @N 5 4 @N 5 (3.28)
@⌘ @⌘ @⌘ @y @y
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]:
Z 1 Z 1 ngp ngp
X X
f (⇠, ⌘) |J| d⇠d⌘ ⇠
= Wi Wj f (⇠i , ⌘j ) |J| . (3.29)
1 1 i=1 j=1
Wi are the weight factors, ⇠i and ⌘i are the Gauss points and ngp is the number
of Gauss integration points, which specify that a polynomial of degree 2ngp 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 Z
T
([B] )⇢cu [N ] dV (3.30)
V
30
Convection modelling
in the system matrix [K], where a three-point Gauss-Lobatto quadrature rule is ap-
plied. 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 bound-
ary 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]
Type of Gauss quadrature Gauss points Weight factors
p
Gauss-Legendre, ngp = 2 ±1/ 3 1
Gauss-Lobatto, ngp = 3 0 4/3

±1 1/3
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 solu-
tions. 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 doc-
umentation 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.
31
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

0 2 1
20 1 0 13 0 1 @u @uv
B +
u
(u r)u = 4@ A @
@/@x u
A5 @ A = B @x @y C
C (3.31)
@ @uv @v 2 A
v @/@y v +
@x @y
and the term r (r u) can be written as
0 1
@u @u
+
B 2 @y 2 C .
r (r u) = u = @ @x
@v @v A (3.32)
+
@x2 @y 2
The Navier-Stokes equation (3.4) and continuity equation in two space dimensions
thus take the form
✓ ◆
@u 1 @p @u2 @uv µ @u @u fx
+ = + 2
+ 2
+ and (3.33)
@t ⇢0 @x @x @y ⇢0 @x @y ⇢0
✓ ◆
@v 1 @p @v 2 @uv µ @v @v fy
+ = + + + , (3.34)
@t ⇢0 @y @y @x ⇢0 @x2 @y 2 ⇢0
@u @v @⇢
+ = =0, (3.35)
@x @y @t
whereby fx and fy 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 signi-
ficantly simpler to work with. The solution is then found by following a three step
approach, which is outlined in the following, where ut denotes the solution at the time
t and ut+1 the final solution at the time t + 1, while u⇤ and u⇤⇤ denote intermediate
solutions of the respective split equations.
32
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.
✓ 2 ◆
@ut @ut vt
u⇤ u t t + =0 (3.36)
@x @y
✓ ◆
⇤ @ut vt @u2t
v vt t + =0 (3.37)
@x @y
This, however, introduces a Courant-Friedrichs-Lewy (CFL) condition, which
limits the time step by a constant times the spatial resolution [38]. For a first-
order upwind scheme, this CFL condition requires that, for stability, the condi-
tions
|u|max t
CF L(x) ⌘ 1 (3.38)
x
and
|v|max t
CF L(y) ⌘ 1 (3.39)
y
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⇤⇤ fx
u⇤⇤ u⇤ = t 2
+ 2
+ t and (3.40)
⇢0 @x @y ⇢0
✓ ◆
⇤⇤ ⇤ µ @v⇤⇤ @v⇤⇤ fy
v v = t + + t . (3.41)
⇢0 @x2 @y 2 ⇢0
3. Pressure correction step

The intermediate velocity field (u⇤⇤ , v⇤⇤ ) is corrected by the gradient of pressure
pt+1 to ensure incompressibility.
t @pt+1
ut+1 u⇤⇤ = (3.42)
⇢0 @x
t @pt+1
vt+1 v ⇤⇤ = (3.43)
⇢0 @y
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
rut+1 = 0 leads to:
⇢0
pt+1 = ru⇤⇤ (3.44)
t
The pressure correction step can therefore be implemented as follows:
(a) Compute the gradient ru⇤⇤ .

⇢0 ⇤⇤
(b) Solve the Poisson equation pt+1 = t ru for pt+1 .
(c) Compute the pressure gradient rpt+1 .

t
(d) Update the velocity field ut+1 = u⇤⇤ ⇢0 rpt+1 .
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
34
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,
@u
as illustrated here for ui,j and at the grid point (i, j):
@x i,j
✓ ◆ ✓ ◆
@2u @2u ui 1,j 2ui,j + ui+1,j ui,j 1 2ui,j + ui,j+1
ui,j = + ⇡ +
@x2 i,j @y 2 i,j
2
x
2
y
✓ ◆ (3.45)
@u ui+1,j ui 1,j
⇡ . (3.46)
@x i,j 2 x
This approximation for the first derivative can yield instabilities [36], but the centered
stencil ✓ ◆
@u ui+1,j ui,j
⇡ (3.47)
@x i+ 12 ,j x
@u
is a stable approximation for between two grid points [36]. This aids the staggered
@x
grid method, where this position is exactly the position of the pressure pi,j .
35
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
✓ ✓ ◆ ◆
= min 1.2 t max max |ui,j | , max |vi,j | , 1 . (3.48)
i,j i,j
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
ui+1,j + ui,j
ūhi+ 1 ,j = and (3.49)
2 2
ui,j+1 + ui,j
ūvi,j+ 1 = , (3.50)
2 2
(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̃h , ũv , etc., the nonlinear terms in equations (3.36) and (3.37) can be written as the
linear combinations
✓ 2 ◆
@u @uv @ ⇣ h 2 ⌘ @
+ = ū ūh ũh + ūv v̄ h v̄ h ũv and (3.51)
@x @y @x @y
✓ ◆
@uv @u2 @ @ ⇣ v 2 ⌘
+ = ūv v̄ h |ūv | ṽ h + (v̄ ) |v̄ v | ṽ v , (3.52)
@x @y @x @y
which are again approximated by centered finite di↵erence stencils.
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 con-
ditions 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
ui,j = 0, if i = 1 and nx + 1 (3.53)
36
at the west and east domain borders and
vi,j = 0, if j = 1 and ny + 1 (3.54)
at the south and north borders. To force zero velocity at the boundary, the conditions
ui,j=1 = ui,j=2 and

(3.55)
ui,j=ny +1 = ui,j=ny +2
at the south and north domain border and
vi=1,j = vi=2,j and
(3.56)
vi=nx +1,j = vi=nx +2,j
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:
pi=1,j = pi=2,j
pi=nx +1,j = pi=nx +2,j

(3.57)
pi,j=1 = pi,j=2
pi,j=ny +1 = pi,j=ny +2 .
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]:
(1 a)2
f= 108 u⇤⇤ (3.58)
a3 + 0.001
37
The large number 108 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 im-
plemented in the Laplace-operator to solve the Poisson equation (3.44). There-
fore, the Laplace-operator must be assembled anew in every time step.
38
Chapter 4
MATLAB implementation
This chapter deals with the final form of the MATLAB implementation of the nu-
merical 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 difficult 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
39
Chapter 4. MATLAB implementation
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 paramet-
ers 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.
40
Structure of the program
START: Run PCM_solver with input datafile
Preprocessing
) Load input data
) Call mesh2d(param) to generate mesh
) Initialize data arrays (T_old, U_old, V_old, ..)
) Call shape_mat to calculate shape matrices
) Preprocess FEM data arrays by calling FEM_preprocessor
Call heat2Delem_new to compute FEM matrices [C], [K] and [R]
Call solve_euler_back to solve FEM equation system and obtain T_new
NO Y ES
Model with convection?
Preprocess velocity and staggered grid data
Call solve_navierstokes to solve Navier-Stokes

equation and obtain U_new, V_new
Postprocess velocity data
Write to data saving arrays
Initialize data arrays for next time step

) T_old = T_new
) U_old = U_new
) V_old = V_new
NO Y ES
Max. number of time steps reached?
Postprocessing
) Save results to output file
) Display information for user
END
Figure 4.1: Flowchart diagram of FEM model PCMsolver

Note. Code given in A.1
41
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 ) ;
42
Description of the key functions
in PCMsolver. As already mentioned in Section 3.2.1, the two-dimensional domain is

split into an equidistant mesh of nel = nx ny four-noded bilinear rectangular finite
elements, which yields the total sum of nn = (nx + 1) (ny + 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 nx = 2 and ny = 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:
[ g a u s s ] = shape mat ( param ) ;
43
[ g a u s s ] = 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 as-
sembly. It evaluates the values of the density, specific heat capacity and heat con-
ductivity for every Gauss point of every element in the domain and preprocesses
boundary informations, as already explained in 4.1. In addition, sparse matrix as-
sembly relations are defined, to be used later in heat2Delem new. The sparse matrix
format in MATLAB provides efficient 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 im-
plementation 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 e l , 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
44
Description of the key functions
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 o l v e e u l e r b a c k ( T old , C, K, R, param ) ;
It contains only one definition:
Tnew = Told + param . mesh . d t ...

( (C + param . mesh . d t K) \ (R K 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 efficient 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 o l v e n a v i e r s t o k e s ( U old , V old , T i j , p e l , param

);
Therein, the force terms, which inhibit fluid motion in the solid elements, are pre-
processed, 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].
45
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 ( i n p u t d a t a f i l e ) ;
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 recom-
mended. 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 re-
spective 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
46
Usage information
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.5 e -4; % thermal expansion coeff . [1 / K ]
74 param . pcm . grav = 9.81; % standard gravity [ m / s ^2]
75 param . pcm . mue = 5.8 e -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 )
47
101 param . air . Tm = 1000; % melting temperature [ C]

102 param . air . dTm = 0; % area of melting [ C]
103 param . air . c_s = 1000; % heat capacity of solid [ J / kg * K ]
104 param . air . lh = 0; % latent heat of PCM [ J / kg ]
105 param . air . k_s = 0.024; % heat conductivity , solid [ W / m * K ]
106
107 %% data recovery
108 param . save . intervall = 6* ceil (1/ param . mesh . d_t ) ;
% time step intervall
109 % for T ,U , V to be saved
110 param . save . make_snapshot = 1; % (0 = off , 1 = on ) % dump snapshots of data
111 param . save . snapshot_int = 120; % time step intervall to make snapshots
112
113 %% display information for user
114 % (0 = off , 1 = on )
115 param . display . waitbar = 1; % show waitbar during run
116 % info : each waitbar update needs about 0.2 seconds computation time
117 param . display . update = 10; % minimum time between waitbar updates in seconds
118 param . display . Tplot = 0; % plot temperature distribution at the end
119 param . display . Vplot = 0; % plot velocity distribution at the end
120 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 under-
stood 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
48
Usage information
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 do-
main. 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 spe-
cifiers 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 coefficient 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.
49
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 num-
ber 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 MATLAB-
internal 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.
50
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 per-
formed. 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 scen-
ario 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.
51
Chapter 5. Validation
5.1 Validation of heat conduction model
To validate the accuracy and convergence of the model without convection, the one-
dimensional Stefan problem was considered, as described in Section 2.1.4. The mater-
ial parameters used in this simulation are those of KNO3 -NaNO3 given in Table 3.1,
but to simplify the analytical solution, the specific heat capacity cl = cs and thermal
conductivity kl = ks were set to the values of the solid substance. The melting point
was defined as Tm = 220 C and the left wall input temperature as 235 C. The initial
temperature of the simulated domain was naturally set to Tinit = Tm ✏, 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]

v
u n
u1 X
RM SE = t
2
(sanalytical ssimulations ) (5.1)
n i=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.
52
Validation of heat conduction model
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)
x
Mesh size 5 2 1 0.5 0.25 0.1
mm
Number of elements 10 25 50 100 200 500
RMSE
(analytical) 9.36 2.03 1.86 2.01 2.02 2.03
10 4 mm
RMSE
(convergence) 11.0 2.07 0.76 0.07 0.02 -
10 4 mm
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.
53
Table 5.2: RMSE of melting front position for di↵erent values of time step size com-
pared to Stefan solution (analytical) and t = 0.1 s simulation (convergence)
t
Time step size 10 5 1 0.5 0.1
s
RMSE
(analytical) 3.66 2.16 2.01 2.10 2.12
10 4 mm
RMSE
(convergence) 4.34 2.19 0.33 0.08 -
10 4 mm
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)
✏
Mushy region size 4 2 1 0.5 0.2 0.1
C
RMSE
(analytical) 8.00 4.04 2.01 1.06 3.04 8.28
10 4 mm
54
Validation of heat conduction model
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.
0.035
0.03
Position of the melting front [m]
0.025
0.02
0.015
analytical solution
0.01 = 4.0°C
= 2.0°C
= 1.0°C
0.005 = 0.5°C
= 0.2°C
= 0.1°C
0
0 1 2 3 4 5
time [h]
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
55
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.
4
The RMSE value compared to the analytical solution was found to be 1.08 10 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.
56
Validation of convection model
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 = 1010 W/mK.
The gallium used in the experiments [5] had a purity of 99.6 % and a melting temper-
ature 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
57
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
v
u n
u1 X
RM SE = t
2
(aexperiment asimulations ) (5.2)
n i=1
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.
58
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 (conver-
gence)
x
Mesh size 4 2 1.6 1 0.8 0.5
mm
Number of elements 352 1408 2200 5632 8800 22528
RMSE
(experimental) 4.72 2.32 2.26 2.31 2.16 1.95
100
RMSE
(convergence) 2.84 0.53 0.39 0.38 0.25 -
100
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)
t
Time step size 50 25 15 10 5 2
ms
RMSE
(experimental) 2.62 2.31 1.87 1.50 1.12 1.34
100
RMSE
(convergence) 2.77 2.63 2.22 1.65 0.94 -
100
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 experi-
mental data. Still, a root mean square error of below 2 %, which accounts to a relative
59
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 com-
pared to test case (experimental)
✏
Mushy region size 0.5 0.25 0.1 0.05 0.02
C
RMSE
(experimental) 1.58 1.14 2.31 2.52 9.70
100
E↵ect of the aspect ratio
Another question with respect to stability and convergence of the developed model
x y
was whether large aspect ratios y or x 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
60
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.
70
legend, including RMSE
experimental data
60
ar: 1/1 : 1.95 %
ar: 1.6/1 : 2.23 %
ar: 2/1 : 2.38 %
melted volume fraction [%]
50 ar: 3.2/1 : 2.60 %

ar: 4/1 : 2.74 %
ar: 8/1 : 4.00 %
40 ar: 8/8 : 4.72 %
30
20
10
0
0 2 4 6 8 10 12 14 16 18 20
time [min]
x
Figure 5.2: Time evolution of liquid fraction for di↵erent aspect ratio values ar =
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
61
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 sim-
ulations 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.
70
legend, including RMSE
experimental data
60
ar: 1/1 : 1.95 %
ar: 1/1.6 : 2.15 %
ar: 1/2 : 1.84 %
50 ar: 1/3.2 : 2.49 %

ar: 1/4 : 3.33 %
ar: 1/8 : 4.82 %
40 ar: 8/8 : 4.72 %
30
20
10
0
0 2 4 6 8 10 12 14 16 18 20
time [min]
x
Figure 5.3: Time evolution of liquid fraction for di↵erent aspect ratio values ar =
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
62
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.
0.06
0.05
0.04
y [m]
0.03
experimental data
simulation t = 5 ms
0.02
simulation t = 50 ms
3 min
6 min
0.01 10 min
15 min
19 min
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x [m]
Figure 5.4: Melting fronts of two simulations with di↵erent time step size compared
to experimental data, given at five time points
63
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.
0.06
0.05
0.04
y [m]
0.03
experimental data
simulation x = 0.5 mm
0.02
simulation x = 2 mm
3 min
6 min
0.01 10 min
15 min
19 min
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x [m]
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 discrep-
ancy 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
64
Cross-validation with ANSYS Fluent
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.
65
(a) after 15 minutes
(b) after 30 minutes
(c) after 60 minutes
(d) after 90 minutes
(e) after 120 minutes
(f) after 150 minutes
(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).
66
Test of a fluctuating charging/discharging scenario
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 gal-
lium 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 man-
ner, 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 max-
imum 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 nu-
merical difficulties 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
67
convective heat transfer leads to the separated zones merging very quickly and there-
fore 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 melt-
ing 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.
50 0.035
45 0.03
0.025
40
max. velocity [m/s]

temperature [°C]
0.02
35
0.015
30
0.01
25
input temperature
mean temperature of gallium cavity 0.005
melting temperature
20 max. velocity of molten gallium
0
0 5 10 15 20 25 30 35 40
time [min]
Figure 5.7: Temperature and velocity evolution over time in fluctuating gallium test
case charging/discharging scenario
68
Test of a fluctuating charging/discharging scenario
0.06
46
0.05 44
42
0.04
temperature [°C]
40
y [mm]
0.03 38
36
0.02
34
0.01
32
30
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x [mm]
Figure 5.8: Temperature distribution in the gallium cavity at the simulation time of
26 minutes during charging/discharging scenario
m/s
0.06
0.025
0.05
0.02
0.04
y [mm]
0.015
0.03
0.01
0.02
0.005
0.01
0
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x [mm]
Figure 5.9: Velocity distribution in the gallium cavity at the simulation time of 26
minutes during charging/discharging scenario
69
70
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 con-
ducted 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. There-
fore, 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 aluminium-
framed 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
71
Chapter 6. Parameter studies
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.
Yaluminium
PCM
Aluminum
Steel
YPCM
Yaluminium
Xsteel XPCM
Xaluminium
Figure 6.1: Sketch of a oftypical

Vienna University TechnologyPCM Dominik
cavity fin section attached to the
Pernsteiner 1 steel contain-
ment 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 mW
2 K . The heat transfer through the well insulated outwards facing wall
was described by the thermal conductivity ↵out = 0.02 mW

2 K . Again, these paramet-
ers 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
72
Study of di↵erent cavity designs
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 geo-
metric aspects of a fin section, meaning the dimensions of the PCM XP CM and YP CM
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 XP CM , YP CM = 5 mm and XP CM , YP CM = 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
YP CM
mm
50 25 10 5
50 10v10 10v20 10v30 10v40
XP CM 25 20v10 20v20 20v30 20v40

mm
10 30v10 30v20 30v30 30v40
5 40v10 40v20 40v30 40v40
73
Table 6.2: Simulated parameter study cases and associated tags with fin width
Yaluminium = 2.5 mm
YP CM
mm
50 25 10 5
50 11v10 11v20 11v30 11v40
XP CM 25 21v10 21v20 21v30 21v40

mm
10 31v10 31v20 31v30 31v40
5 41v10 41v20 41v30 41v40
6.1.2 Simulation results
To analyse the simulation results and compare the di↵erent cavity dimensions, the
H
evolution of the enthalpy change per surface area was calculated, where A is
A
the surface area between the PCM fin section and the steam storage. The surface
therefore takes the value
A = (YP CM + 2Yaluminium ) 1m . (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 (XP CM = 50 mm, 25 mm),
the lower cavities can store more enthalpy per surface area and the charging and dis-
charging 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
74
106 106
8 8
total enthalpy total enthalpy
PCM latent heat PCM latent heat
7 7
enthalpy per surface area [J/m2 ]

case 10v10 case 20v10
6 6
5 5
4 4
3 3
2 2
1 1
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
time [h] time [h]
(a) cases with constant XP CM = 50 mm (b) cases with constant XP CM = 25 mm
106 106
8 8
7 7

6 6
5 5
4 4
3 3
2 2
1 1
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
time [h] time [h]
(c) cases with constant XP CM = 10 mm (d) cases with constant XP CM = 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 di-
mensioning of the cavity is also a question of cost. If unused PCM material must be
75
106 106
8 8
7 7

6 6
5 5
4 4
3 3
2 2
1 1
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
time [h] time [h]
(a) cases with constant XP CM = 50 mm (b) cases with constant XP CM = 25 mm
106 106
8 8
7 7

6 6
5 5
4 4
3 3
2 2
1 1
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
time [h] time [h]
(c) cases with constant XP CM = 10 mm (d) cases with constant XP CM = 5 mm
Figure 6.3: Enthalpy per surface area of the simulated cases with fin width
Yaluminium = 2.5 mm
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.
76
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 dimin-
ished 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.
106 106
8 8
7 7
PCM latent heat enthalpy per surface area [J/m2 ] PCM latent heat
6 6
5 5
4 4
3 3
2 2
1 1
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
time [h] time [h]
(a) cases with constant XP CM = 50 mm and (b) cases with constant XP CM = 5 mm and
YP CM = 5 mm YP CM = 50 mm
Figure 6.4: Enthalpy per surface area of the simulated cases with di↵erent fin width
Yaluminium
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
77
in Section 6.2.2 with XP CM = 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 cross-
validation 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
78
Comparison of di↵erent cavity orientations
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 ca-
pacity 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 XP CM = 118 mm and YP CM = 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 mW
2 K . 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 temper-
ature 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.
79
HyStEPs: Location of a Fin Section
0°
𝐺Ԧ
Steam
90°
Storage
180°
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
80
Comparison of di↵erent cavity orientations
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.
105
7
5
enthalpy [J]
3
total enthalpy
PCM latent heat
2 0 degrees
45 degrees
90 degrees
1 135 degrees
180 degrees
conduction only
0
0 0.5 1 1.5 2 2.5 3
time [h]
Figure 6.6: Change of enthalpy during the charging process

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.
81
100 100
0 degrees
90 90 45 degrees
90 degrees
80 80
135 degrees

180 degrees
70 70
conduction only
60 60
50 50
40 40
0 degrees
30 45 degrees 30
90 degrees
20 20
135 degrees
180 degrees
10 10
conduction only
0 0
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
time [h] time [h]
(a) charging (b) discharging
Figure 6.7: Evolution of the melted volume fraction during charging (a) and dischar-
ging (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.
82
Chapter 7
Conclusion
This last chapter presents a short summary of the overall approach to the stated prob-
lems 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 con-
vection. 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,
83
Chapter 7. Conclusion
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 mater-
ials 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 pre-
dictions 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.
84
Scope and limitations
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 con-
vection 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 para-
meter 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 difficulties 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 para-
meters. 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
85
change and would not be exposed to the e↵ect of the latent heat. Even if these limit-
ations 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 coefficient. Such an approach was followed for some time to further im-
prove 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 ap-
proach 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 con-
trol 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 mater-
ial 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 computa-
tional 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 sufficient approxim-
ation of the discharging process by only considering heat conduction was also reported
in reference [27] for similar PCM cavity geometry.
86
Suggested future objectives
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 un-
derway to provide a detailed foundation for design optimization of the considered
PCM cell geometry [49]. Therein, varying aluminium proportions, varying fin spa-
cings 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.
87
Furthermore, adaptive mesh or even moving mesh discretization methods, as docu-

mented 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 at-

tempted 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 explan-
ation 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
@(⇢cT )
= r(krT ) r(⇢cT u) + S . (7.1)
@t
In this form however, c takes a fixed value instead of the apparent heat capacity, and
the source term
a
S= H (7.2)
t
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
n o
[C] Ṫ + [K] {T } = [R] . (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:
88
Suggested future objectives
1. The source term (7.2) is approximated by using an estimate of the liquid fraction
field a.
2. Equation (3.23) can be evaluated and solved by a backwards Euler step.
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= Tm
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].
89
90
Bibliography
[1] R. Hofmann, S. Dusek, S. Gruber and G. Drexler-Schmid. Design Optimization
of a Hybrid Steam-PCM Thermal Energy Storage for Industrial Applications.
Energies, 12(5), 2019. doi:10.3390/en12050898.
[2] S. Gruber. Charakterisierung eines Hybridspeichers zur Dampfbereitstellung

in industriellen Prozessen. Theoretische Untersuchung und Konzeption eines
Pruefstandes. Master’s thesis, TU Wien, Vienna, Austria, 2018. URL http:
//katalog.ub.tuwien.ac.at/AC15194551.
[3] Y. Dutil, D.R. Rousse, N. Salah, S. Lassue and L. Zalewski. A review on phase-
change materials: Mathematical modeling and simulations. Renewable and Sus-
tainable Energy Reviews, 15:112–130, 2011. doi:10.1016/j.rser.2010.06.011.
[4] MATLAB. version 9.3.0 (R2017b). The MathWorks Inc., Natick, Massachusetts,
2017.
[5] C. Gau and R. Viskanta. Melting and Solidification of a Pure Metal on a Vertical
Wall. Journal of Heat Transfer, 108, 02 1986. doi:10.1115/1.3246884.
[6] A. Kumar and S.K. Shukla. A Review on Thermal Energy Storage Unit for Solar
Thermal Power Plant Application. Energy Procedia, 74:462 – 469, 2015. ISSN
1876-6102. doi:10.1016/j.egypro.2015.07.728.
[7] K.K. Pillai and B.J. Brinkworth. The storage of low grade thermal energy using
phase change materials. Applied Energy, 2(3):205 – 216, 1976. ISSN 0306-2619.
doi:10.1016/0306-2619(76)90025-8.
[8] P. W. Egolf and H. Manz. Theory and modeling of phase change materials with
and without mushy regions. International Journal of Heat and Mass Transfer,
37(18):2917 – 2924, 1994. ISSN 0017-9310. doi:10.1016/0017-9310(94)90346-8.
91
Bibliography
[9] B. Nedjar. An enthalpy-based finite element method for nonlinear heat problems
involving phase change. Computers and Structures, 80(1):9 – 21, 2002. ISSN
0045-7949. doi:10.1016/S0045-7949(01)00165-1.
[10] W. Demtröder. Experimentalphysik 1: Mechanik und Wärme. Springer-

Lehrbuch. Springer Berlin Heidelberg, 2018. doi:10.1007/978-3-662-54847-9.
[11] R.T. Tenchev, J.A. Mackenzie, T.J. Scanlon and M.T. Stickland. Finite ele-
ment moving mesh analysis of phase change problems with natural convec-
tion. International Journal of Heat and Fluid Flow, 26(4):597 – 612, 2005.
doi:10.1016/j.ijheatfluidflow.2005.03.003.
[12] S. Chandrasekhar. Stochastic problems in physics and astronomy. Reviews of

modern physics, 15:1–89, 1943. doi:10.1103/RevModPhys.15.1.
[13] S. Braun. Strömungslehre für TPh. TU Wien, Institut für Strömungsmechanik

und Wärmeübertragung, 2017.
[14] P. Dawkins. The Heat Equation. Available at http://tutorial.math.lamar.edu/

Classes/DE/TheHeatEquation.aspx.
[15] T. Jonsson. On the one dimensional stefan problem : with some numerical
analysis, 2013. URL http://www.diva-portal.org/smash/record.jsf?pid=diva2%
3A647481&dswid=-7335.
[16] T.F. Cheng. Numerical analysis of nonlinear multiphase Stefan problems. Com-
puters & Structures, 75(2):225 – 233, 2000. ISSN 0045-7949. doi:10.1016/S0045-
7949(99)00071-1.
[17] S. Liu, Y. Li and Y. Zhang. Mathematical solutions and numerical mod-

els employed for the investigations of PCMs phase transformations. Renew-
able and Sustainable Energy Reviews, 33:659 – 674, 2014. ISSN 1364-0321.
doi:10.1016/j.rser.2014.02.032.
[18] H.J. Stetter. Analysis of Discretization Methods for Ordinary Di↵erential Equa-
tions. Springer Tracts in Natural Philosophy. Springer Berlin Heidelberg, 2013.
doi:10.1002/zamm.19750550213.
92
Bibliography
[19] D.W. Pepper and J.C. Heinrich. The Finite Element Method: Basic Concepts
and Applications. Series in Computational and Physical Processes in Mechanics
and Thermal Sciences. CRC Press, 2005. doi:10.1201/9780203942352.
[20] M. Kaltenbacher. Numerical Simulation of Mechatronic Sensors and Actuators:

Finite Elements for Computational Multiphysics, third edition. Springer Berlin
Heidelberg, 01 2015. doi:10.1007/978-3-642-40170-1.
[21] Mathematik Uni Dortmund. Finite di↵erence method. Available at http://www.

mathematik.uni-dortmund.de/˜kuzmin/cfdintro/lecture4.pdfx.
[22] P.D.M. Causon, P.C.G. Mingham and L. Qian. Introductory Finite Volume Meth-
ods for PDEs. Bookboon, 2011. ISBN 978-87-7681-882-1.
[23] F. Moukalled, L. Mangani and M. Darwish. The Finite Volume Method in Com-
putational Fluid Dynamics: An Advanced Introduction with OpenFOAM and
Matlab. Fluid Mechanics and Its Applications. Springer International Publishing,
2015. doi:10.1007/978-3-319-16874-6.
[24] R.D. Richtmyer and K.W. Morton. Di↵erence methods for initial-value problems.
Interscience tracts in pure and applied Mathematics. Interscience Publishers,
1967. doi:10.1137/1010073.
[25] S. Dusek and R. Hofmann. A Hybrid Energy Storage Concept for Future Ap-
plication in Industrial Processes. Thermal Science, 22(5):2235 – 2242, 2018.
doi:10.2298/TSCI171230270D.
[26] M. M. Farid, A. M. Khudhair, S. A. K. Razack and S. Al-Hallaj. A re-

view on phase change energy storage: materials and applications. Energy
Conversion and Management, 45(9):1597 – 1615, 2004. ISSN 0196-8904.
doi:10.1016/j.enconman.2003.09.015.
[27] J. Vogel, M. Johnson, M. Eck and D. Laing. Numerical analysis of natural

convection in a latent heat thermal energy storage system containing rectangular
enclosures. 2014. URL https://elib.dlr.de/113475/.
93
Bibliography
[28] J. Lopez, Z. Acem and E. P. Del Barrio. KNO3/NaNO3 - Graphite materials for
thermal energy storage at high temperature: Part II. - Phase transition proper-
ties. Applied Thermal Engineering, 30(13):1586 – 1593, 2010. ISSN 1359-4311.
doi:10.1016/j.applthermaleng.2010.03.014.
[29] D.R. Lide. CRC Handbook of Chemistry and Physics, 85th Edition. Number
Bd. 85 in CRC Handbook of Chemistry and Physics, 85th Ed. Taylor & Francis,
2004. ISBN 9780849304859.
[30] S. Dusek and R. Hofmann. Modeling of a Hybrid Steam Storage and Validation
with an Industrial Ruths Steam Storage Line. Energies, 12(6), 2019. ISSN 1996-
1073. doi:10.3390/en12061014.
[31] D. Pernsteiner, L. Kasper, A. Schirrer, R. Hofmann and S. Jakubek.

Co-simulation methodology of a hybrid latent-heat thermal energy storage
unit. Applied Thermal Engineering, page 115495, 2020. ISSN 1359-4311.
doi:10.1016/j.applthermaleng.2020.115495.
[32] H. Herr. Wärmelehre: Technische Physik Band 3. Europa-Lehrmittel, 4 edition,

2006. ISBN 978-3-8085-5064-9.
[33] M.J. Huang, P.C. Eames and B. Norton. Comparison of a small-scale 3D PCM
thermal control model with a validated 2D PCM thermal control model. Solar
Energy Materials and Solar Cells, 90(13):1961 – 1972, 2006. ISSN 0927-0248.
doi:10.1016/j.solmat.2006.02.001.
[34] C. H. Stolze. A history of the divergence theorem. Historia Mathematica, 5(4):

437 – 442, 1978. ISSN 0315-0860. doi:10.1016/0315-0860(78)90212-4.
[35] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions: With

Formulas, Graphs, and Mathematical Tables. Applied mathematics series. Dover
Publications, 1965. ISBN 9780486612720.
[36] B. Seibold. A compact and fast Matlab code solving the incompressible Navier-
Stokes equations on rectangular domains, 2018. Available at http://math.mit.
edu/˜gs/cse/codes/mit18086 navierstokes.pdf.
94
Bibliography
[37] H. P. Langtangen, K. A. Mardal and R. Winther. Numerical methods for in-

compressible viscous flow. Advances in Water Resources, 25:1125–1146, 08 2002.
doi:10.1016/S0309-1708(02)00052-0.
[38] R. Courant, K. Friedrichs and H. Lewy. Über die partiellen Di↵erenzengleichun-

gen der mathematischen Physik. Mathematische Annalen, 100(1):32–74, 1928.
doi:10.1007/BF01448839.
[39] Massachusetts Institute of Technology. Stability of Finite Di↵erence Methods.

Available at http://web.mit.edu/16.90/BackUp/www/pdfs/Chapter14.pdf.
[40] M. Griebel, T. Dornsheifer and T. Neunhoe↵er. Numerical Simulation in Fluid

Dynamics: A Practical Introduction. Monographs on Mathematical Model-
ing and Computation. Society for Industrial and Applied Mathematics, 1998.
doi:10.1137/1.9780898719703.
[41] J. Dallaire and L. Gosselin. Numerical modeling of solid-liquid phase change

in a closed 2D cavity with density change, elastic wall and natural convection.
International Journal of Heat and Mass Transfer, 114:903 – 914, 2017. ISSN
0017-9310. doi:10.1016/j.ijheatmasstransfer.2017.06.104.
[42] GP Nikishkov. Introduction to the Finite Element Method. University of Aizu,

2009. Available at http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.
521.9293&rep=rep1&type=pdf.
[43] K. Atkinson and W Han. Elementary numerical analysis. John Wiley & Sons,
1993. ISBN 9780471897330.
[44] A. D. Brent, V. R. Voller and K. J. Reid. Enthalpy-porosity technique for model-

ing convection-di↵usion phase change: Application to the melting of a pure metal.
Numerical Heat Transfer, 13(3):297–318, 1988. doi:10.1080/10407788808913615.
[45] M. Breuer. A challenging test case for large eddy simulation: high Reynolds
number circular cylinder flow. International Journal of Heat and Fluid Flow, 21
(5):648 – 654, 2000. ISSN 0142-727X. doi:10.1016/S0142-727X(00)00056-4.
95
Bibliography
[46] Z. Boz, F. Erdogdu and M. Tutar. E↵ects of mesh refinement, time step size
and numerical scheme on the computational modeling of temperature evolution
during natural-convection heating. Journal of Food Engineering, 123:8 – 16, 2014.
ISSN 0260-8774. doi:10.1016/j.jfoodeng.2013.09.008.
[47] J. Vogel, J. Felbinger and M. Johnson. Natural convection in high temperature

flat plate latent heat thermal energy storage systems. Applied Energy, 184:184–
196, 2016. doi:10.1016/j.apenergy.2016.10.001.
[48] W.B. Ye. Thermal and hydraulic performance of a rectangular stor-

age cavity. Applied Thermal Engineering, 93:1114 – 1123, 2016.
doi:10.1016/j.applthermaleng.2015.10.083.
[49] L. Kasper, D. Pernsteiner, M. Koller, A. Schirrer, S. Jakubek and R. Hofmann.

Numerical studies on the melting and solidification process in an rectangular
aluminium finned latent-heat thermal energy storage. in submission. Applied
Thermal Engineering, 2020. Manuscript submitted for publication, under review.
[50] D. Kavetski, P. Binning and S. W. Sloan. Adaptive backward Euler time stepping
with truncation error control for numerical modelling of unsaturated fluid flow.
International Journal for Numerical Methods in Engineering, 53(6):1301–1322,
2002. doi:10.1002/nme.329.
[51] V.R. Voller and C.R. Swaminathan. Eral Source-Based Method for Solidification
Phase Change. Numerical Heat Transfer Part B - Fundamentals, 19, 04 1991.
doi:10.1080/10407799108944962.
[52] V.R. Voller, C.R. Swaminathan and B. Thomas. Fixed Grid Techniques for Phase
Change Problems: A Review. International Journal for Numerical Methods in
Engineering, 30:875 – 898, 09 1990. doi:10.1002/nme.1620300419.
96
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
97
Appendix A. MATLAB code listings
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
liq uid ific atio n
69 Tl = param . pcm . Tm - param . pcm . dTm ; % temperature of partial
liq uid ific atio n
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 ] = FE M _p r ep r oc es s or ( 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 = s o lv e _e u le r _b a ck ( 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
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 ] = s o l v e _ n a v i e r s t o k e s ( 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) ) ;
98
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
155 [f , timearray , up ] = wa i tb a r2 D _s i mp l e (f , param , timearray , toc ,t , up ) ;
156 end
157 %% dump variables
158 if mod (t , param . save . intervall ) == 0
159 Tt (: , s ) = T_new ;
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 %03 d . mat , param . save . name
,t);
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
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
99
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 ;
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
100
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
19 gauss . B2 = BmatHeat2D ( gp (2) , gp (1) , coord ) ;
21
26
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 ) ;
59
60 gauss . Nc3 = NmatHeat2D ( gp (1) , gp (3) ) ; % eta -1 , psi 1
64
65 gauss . Nc4 = NmatHeat2D ( gp (2) , gp (1) ) ; % eta 0 , psi -1
101

69
70 gauss . Nc5 = NmatHeat2D ( gp (2) , gp (2) ) ; % eta 0 , psi 0
74
79
80 gauss . Nc7 = NmatHeat2D ( gp (3) , gp (1) ) ; % eta 1 , psi -1
84
89
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 ] = FE M _p r ep r oc e ss o r ( 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 ;
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 ;
102
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 (: ,:) ) ;
38 end
39 % gauss point 3
40 T3 = gauss . N3 * T_old ( IEN (: ,:) ) ;
48 end
49 % gauss point 4
50 T4 = gauss . N4 * T_old ( IEN (: ,:) ) ;
58 end
59
60 % density of element
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
68 str1 = strcat ( param . , param . geom . material {m ,2} , . rho ) ;
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;
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.5 7735 027) /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 (: ,:) ) ;
103
95 yl1 = yl (1 ,:) +( yl (4 ,:) - yl (1 ,:) ) *(1 -0.57735027) /2;

96 yl2 = yl (1 ,:) +( yl (4 ,:) - yl (1 ,:) ) *(1 +0.5 7735 027) /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
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 ;
19 hl_1 = 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) ;
104
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 (: ,:) ) ;
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 ) ;
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 ) ;
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) +...
105
114 ( solid ( circshift ( flag_el ==0 , nx ) ) & ...

115 solid ( circshift ( flag_el ==0 , - nx ) ) ) *2;
116 % left side element
117 correct_u ( flag_el ==3) = -1+ solid ( circshift ( flag_el ==3 ,1) ) *3;
123 correct_v ( flag_el ==32) = 1+ solid ( circshift ( flag_el ==32 , - nx ) ) *1;
125 correct_v ( flag_el ==34) = -1+ solid ( circshift ( flag_el ==32 , nx ) ) *3;
126 % right side element
127 correct_u ( flag_el ==1) = 1+ solid ( circshift ( flag_el ==1 , -1) ) *1;
136 % upper side element
142 % lower side element
148
149 % perform velocity correction
150 u =[ u1 , u2 , u3 , u4 , u5 , u6 , u7 , u8 , u9 ];
151 v =[ v1 , v2 , v3 , v4 , v5 , v6 , v7 , v8 , v9 ];
152 u ( correct_u ==2 ,:) = 0;
153 u ( correct_u ==1 ,:) = max (0 , u ( correct_u ==1 ,:) ) ;
154 u ( correct_u == -1 ,:) = min (0 , u ( correct_u == -1 ,:) ) ;
155 v ( correct_v ==2 ,:) = 0;
156 v ( correct_v ==1 ,:) = max (0 , v ( correct_v ==1 ,:) ) ;
157 v ( correct_v == -1 ,:) = min (0 , v ( correct_v == -1 ,:) ) ;
158 u1 = u (: ,1) ; u2 = u (: ,2) ; u3 = u (: ,3) ; u4 = u (: ,4) ; u5 = u (: ,5) ;
159 u6 = u (: ,6) ; u7 = u (: ,7) ; u8 = u (: ,8) ; u9 = u (: ,9) ;
160 v1 = v (: ,1) ; v2 = v (: ,2) ; v3 = v (: ,3) ; v4 = v (: ,4) ; v5 = v (: ,5) ;
161 v6 = v (: ,6) ; v7 = v (: ,7) ; v8 = v (: ,8) ; v9 = v (: ,9) ;
162 end
163
164 %% calculate element matrices
165 Ce2 = gauss . detJ *(...
166 bsxfun (@ times , gauss . N1_q , permute ( gauss . rho . * c1 ,[3 1 2]) ) +...
169 bsxfun (@ times , gauss . N4_q , permute ( gauss . rho . * c4 ,[3 1 2]) ) ) ;
170 C_e = reshape ( Ce2 ,4*4* nel ,1) ;
171
172 KCe = gauss . detJ *(...
173 bsxfun (@ times , gauss . B1_q , permute ( k1 ,[3 1 2]) ) +...
176 bsxfun (@ times , gauss . B4_q , permute ( k4 ,[3 1 2]) ) ) ;
177
178 right = param . mesh . flag (: ,1) == 1; % right side of PCM
179 KB1e = gauss . dS_e *( ...
180 bsxfun (@ times , gauss . Ne1_q , permute ( right . * gauss . alpha_out1 ( t ) ,[3 1
2]) ) +...
181 bsxfun (@ times , gauss . Ne2_q , permute ( right . * gauss . alpha_out2 ( t ) ,[3 1
2]) ) ) ;
182
183 RB1e = gauss . dS_e *( ...
184 bsxfun (@ times , gauss . Ne1 , permute ( right . * gauss . alpha_out1 ( t ) .* gauss .
Tout1 ( t ) ,[3 1 2]) ) +...
185 bsxfun (@ times , gauss . Ne2 , permute ( right . * gauss . alpha_out2 ( t ) .* gauss .
Tout2 ( t ) ,[3 1 2]) ) ) ;
186
187 left = param . mesh . flag (: ,1) == 3; % left side of PCM
188 KB3e = - gauss . dS_w *( ...
106
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]) ) +
...
...
...
...
...
...
...
...
207 bsxfun (@ times , gauss . NBc1v , permute ( p_el_FEM . * v1 , [ 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) ;
107
248 % Ce = rho ( e ) * gauss . detJ *( c1 ( e ) .* gauss . N1_q +...

249 % c2 ( e ) .* gauss . N2_q + c3 ( e ) .* gauss . N3_q +...
250 % c4 ( e ) .* gauss . N4_q ) ;
251 KCe = gauss . detJ *( k1 ( e ) .* gauss . B1_q + ...
252 k2 ( e ) .* gauss . B2_q + k3 ( e ) .* gauss . B3_q +...
253 k4 ( e ) .* gauss . B4_q ) ;
254
255 % calculate enthaply
256 if param . mesh . elem_mat ( e ) == 1
257 H { e } = calc_h1 ( e ) + calc_h2 ( e ) + calc_h3 ( e ) + calc_h4 ( e ) ;
258 H_l { e } = hl_1 ( e ) + hl_2 ( e ) + hl_3 ( e ) + hl_4 ( e ) ;
259 else
260 H { e } = calc_h ( e ) *( T1 ( e ) + T2 ( e ) + T3 ( e ) + T4 ( e ) ) ;
261 end
262
263 if param . mesh . flag (e ,1) == 1 % right side of PCM
264 nodes = param . mesh . IEN (: , e ) ;
266 y1 = coord (1 ,2) +( coord (4 ,2) - coord (1 ,2) ) *( -0.57735027/2) ;
267 y2 = coord (1 ,2) +( coord (4 ,2) - coord (1 ,2) ) *(0 .577 3502 7/2 ) ;
268
269 KB1e = param . icbc . alpha_out * gauss . dS_e *( ...
270 gauss . Ne1_q + gauss . Ne2_q ) ;
271
272 RB1e = param . icbc . alpha_out * gauss . dS_e *...
273 ( param . icbc . Tout (t , y1 ) * gauss . Ne1 + . . .
274 param . icbc . Tout (t , y2 ) * gauss . Ne2 ) ;
275
276 % calculate heat flow
277 Qout { e } = gauss . dS_e * param . icbc . alpha_out * param . mesh . d_t *(...
278 param . icbc . Tout (t , y1 ) - Te1 ( e ) + ...
279 param . icbc . Tout (t , y2 ) - Te2 ( e ) ) ;
280
281 elseif param . mesh . flag (e ,1) == 3 % left side of PCM
282 nodes = param . mesh . IEN (: , e ) ;
284 y1 = coord (1 ,2) +( coord (4 ,2) - coord (1 ,2) ) *( -0.57735027/2) ;
285 y2 = coord (1 ,2) +( coord (4 ,2) - coord (1 ,2) ) *(0 .577 3502 7/2 ) ;
286
287 KB3e = - gauss . dS_w *( param . icbc . alpha_in (t , y1 ) * gauss . Nw1_q ...
288 + param . icbc . alpha_in (t , y2 ) * gauss . Nw2_q ) ;
289
290 RB3e = - gauss . dS_w *( param . icbc . alpha_in (t , y1 ) *...
291 param . icbc . Tin (t , y1 ) * gauss . Nw1 + . . .
292 param . icbc . alpha_in (t , y2 ) * param . icbc . Tin (t , y2 ) ...
293 * gauss . Nw2 ) ;
294
295 % calculate heat flow
296 Qin { e } = - gauss . dS_w * param . mesh . d_t *(...
297 param . icbc . alpha_in (t , y1 ) *( param . icbc . Tin (t , y1 ) - Tw1 ( e ) ) ...
298 + param . icbc . alpha_in (t , y2 ) *( param . icbc . Tin (t , y2 ) - Tw2 ( e ) ) ) ;
299 end
300
301 % convection term
302 if p_el_FEM ( e ) == 1
303
304 KCe = KCe + u1 ( e ) * gauss . NBc1u + v1 ( e ) * gauss . NBc1v +...
305 u2 ( e ) * gauss . NBc2u + v2 ( e ) * gauss . NBc2v +...
312 u9 ( e ) * gauss . NBc9u + v9 ( e ) * gauss . NBc9v ;
313 end
314 K_e { e } = reshape ( KCe + KB1e + KB3e ,4*4 ,1) ;
315 R_e { e } = RB1e + RB3e ;
316 C_e { e } = reshape ( Ce ,4*4 ,1) ;
317 end
318
319 %% assemble to global matrices and vectors
320
321 % C = sparse ( sp_mat (: ,1) , sp_mat (: ,2) , cell2mat ( C_e ) ,non , non ) ;
322 % K = sparse ( sp_mat (: ,1) , sp_mat (: ,2) , cell2mat ( K_e ) ,non , non ) ;
323 % R = sparse ( sp_vec (: ,1) , sp_vec (: ,2) , cell2mat ( R_e ) ,non ,1) ;
324 %}
325 end
108
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
109
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 ] = so l ve _ eu l er _ ba c k ( 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 ] = w ai t ba r 2D _ si mp l e (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
110
Code Listing A.12: Code of function solve navierstokes

1 %% solve navier stokes
2 function [ U_new , V_new , grad ,U , V ] = s o l v e _ n a v i e r s t o k e s (U ,V , T_ij , p_el , param )
3 % Solves the navier stokes equations via finite difference method
4
5 %% preprocessing
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
111
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;
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
112
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.
113
Appendix B. Additional plots
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(a) temperature after 22 minutes (b) velocity after 22 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(c) temperature after 45 minutes (d) velocity after 45 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(e) temperature after 67 minutes (f) velocity after 67 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(g) temperature after 90 minutes (h) velocity after 90 minutes
Figure B.1: Temperature and velocity distribution at di↵erent simulation times for
cavity orientation 0
114
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(i) temperature after 112 minutes (j) velocity after 112 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(k) temperature after 135 minutes (l) velocity after 135 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(m) temperature after 157 minutes (n) velocity after 157 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(o) temperature after 180 minutes (p) velocity after 180 minutes
115
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(q) temperature after 22 minutes (r) velocity after 22 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(s) temperature after 45 minutes (t) velocity after 45 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(u) temperature after 67 minutes (v) velocity after 67 minutes
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
(w) temperature after 90 minutes (x) velocity after 90 minutes
116
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
117
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
118
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
119
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
120
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
121
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
122
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
temperature [°C]
220 222 224 226 228 230
25
20
y [mm]
10
0
0 20 40 60 80 100 120
x [mm]
123
To increase the efficiency of energy- mentation, considering heat transfer by
intensive industrial processes, thermal both conduction and natural convection.
energy storages can offer new possibi- This successfully validated code can
lities. A novel approach is investigated handle any desired layout of materials
in the project HyStEPs. In this concept, arranged on a rectangular domain.
containers filled with PCM are placed at
the shell surface of a Ruths steam storage, Furthermore, a parameter study of diffe-
to increase storage efficiency. rent dimensions and orientations of the
PCM cavity was conducted. The impact
In this work, a two-dimensional model of natural convection was found to lead
using the finite element method is to significantly varying behaviour of the
developed to simulate the PCM of the studied cavities with different orientation
hybrid storage as designed in the HyStEPs during the charging process, while it
project. The apparent heat capacity was found to be negligible during the
method is applied in a MATLAB imple- discharging process.
ISBN 978-3-85448-036-5
www.tuwien.at/academicpress

Modeling of The Phase Change Material of A Hybrid Storage Using The Finite Element Method

Uploaded by

Copyright:

Available Formats

Modeling of The Phase Change Material of A Hybrid Storage Using The Finite Element Method

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Modeling of The Phase Change Material of A Hybrid Storage Using The Finite Element Method

Uploaded by

Copyright:

Available Formats

Lukas Kasper

MODELING OF THE PHASE

MODELING OF THE PHASE

TU Wien Academic Press 2020

c/o TU Wien Bibliothek

ISBN (print): 978-3-85448-036-5

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

Media proprietor: TU Wien, Karlsplatz 13, 1040 Wien

Proofreading: Jennifer Weck

Univ.Prof. Dr. techn. René Hofmann

Furthermore, I want to acknowledge my colleagues from the Control and Process

Finally, I would like to express my gratitude to my parents and my girlfriend Nicki

List of Figures xvii

List of Tables xxi

List of code Listings xxiii

1.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3 Aim of this work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.4 Methodological approach . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.5 Structure of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1 Basics of thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.3 Discretization methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.1 Model equations and assumptions . . . . . . . . . . . . . . . . . . . . . 22

3.3 Convection modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4.1 Structure of the program . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.2 Description of the key functions . . . . . . . . . . . . . . . . . . . . . . 42

4.3 Usage information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.1 Validation of heat conduction model . . . . . . . . . . . . . . . . . . . 52

5.2 Validation of convection model . . . . . . . . . . . . . . . . . . . . . . 57

5.3 Cross-validation with ANSYS Fluent . . . . . . . . . . . . . . . . . . . 64

5.4 Test of a fluctuating charging/discharging scenario . . . . . . . . . . . 67

6.1 Study of di↵erent cavity designs . . . . . . . . . . . . . . . . . . . . . . 71

6.2 Comparison of di↵erent cavity orientations . . . . . . . . . . . . . . . . 78

6.3 Discussion of conduction vs. convection simulations . . . . . . . . . . . 82

7.2 Scope and limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

7.3 Suggested future objectives . . . . . . . . . . . . . . . . . . . . . . . . 87

A MATLAB code listings 97

B Additional plots 113

3.1 Four-noded bilinear element in natural coordinate system ⇠, ⌘ with

3.2 Two-dimensional domain discretized into a staggered grid . . . . . . . 35

4.1 Flowchart diagram of FEM model PCMsolver . . . . . . . . . . . . . . 41

4.2 Image of a 2x2 mesh, with numbered nodes and elements . . . . . . . 43

5.6 Phase structure in PCM cavity at di↵erent times . . . . . . . . . . . . 66

5.7 Temperature and velocity evolution over time in fluctuating gallium

5.8 Temperature distribution in the gallium cavity at the simulation time

5.9 Velocity distribution in the gallium cavity at the simulation time of 26

6.5 Schematic illustration of the five simulated cavity orientations in rela-

6.6 Change of enthalpy during the charging process . . . . . . . . . . . . . 81

B.1 Temperature and velocity distribution at di↵erent simulation times for

B.2 Temperature and velocity distribution at di↵erent simulation times for

B.3 Temperature and velocity distribution at di↵erent simulation times for

B.4 Temperature and velocity distribution at di↵erent simulation times for

B.5 Temperature and velocity distribution at di↵erent simulation times for

A.1 Code of main function PCMsolver . . . . . . . . . . . . . . . . . . . . . 97

A.2 Code of function mesh2d . . . . . . . . . . . . . . . . . . . . . . . . . . 99

A.3 Code of function shape mat . . . . . . . . . . . . . . . . . . . . . . . . 101

A.4 Code of shape matrix function NmatHeat2D . . . . . . . . . . . . . . . 102

A.5 Code of shape matrix function BmatHeat2D . . . . . . . . . . . . . . . 102

A.6 Code of FEM preprocessing function FEM preprocessor . . . . . . . . . 102