2 D Model

See discussions, stats, and author profiles for this publication at: https://www.researchgate.
net/publication/322499834
State and state of charge estimation for a latent heat storage
Article in Control Engineering Practice · March 2018

DOI: 10.1016/j.conengprac.2017.11.006
CITATIONS READS
22 992
7 authors, including:
Tilman Barz Dominik Seliger

AIT Austrian Institute of Technology TU Wien
92 PUBLICATIONS 837 CITATIONS 1 PUBLICATION 22 CITATIONS
SEE PROFILE SEE PROFILE
Klemens Marx Andreas Sommer

AIT Austrian Institute of Technology Universität Heidelberg
16 PUBLICATIONS 176 CITATIONS 11 PUBLICATIONS 113 CITATIONS
SEE PROFILE SEE PROFILE
Some of the authors of this publication are also working on these related projects:
WHole Battery - Wärmeoptimierte Hochleistungs-Batterie View project
AdaMo Project View project
All content following this page was uploaded by Tilman Barz on 15 January 2018.
The user has requested enhancement of the downloaded file.

State and state of charge estimation
for a latent heat storage
Tilman Barz1 , Dominik Seliger, Klemens Marx

AIT Austrian Institute of Technology GmbH, Center for Energy, Giefinggasse 2, 1210
Vienna, Austria
Andreas Sommer, Sebastian F. Walter, Hans Georg Bock

Interdisciplinary Center for Scientific Computing (IWR), Universität Heidelberg,
Im Neuenheimer Feld 205, 69120 Heidelberg, Germany
Stefan Körkel
OTH Ostbayerische Technische Hochschule Regensburg, Faculty of Computer Science and
Mathematics, Postfach 120327, 93025 Regensburg, Germany
Abstract
A nonlinear state observer is designed for a thermal energy storage with
solid/liquid phase change material (PCM). Using a physical 2D dynamic model,
the observer reconstructs transient spatial temperature fields inside the storage
and estimates the stored energy and the state of charge. The observer has been
successfully tested with a lab-scale latent heat storage with a single pass tube
bundle and the phase change material located in a shell around each tube. It
turns out that the observer robustly tracks the real process data with as few as
four internal PCM temperature sensors.
Keywords: latent heat thermal energy storage (LHTES), state of charge
(SOC), reduced model, Orthogonal Collocation, heat conduction in cylindrical
shell, nonlinear state observer, Kalman filter
1 Corresponding author: tilman.barz@ait.ac.at
Post-print version of the article: Barz, T., Seliger, D., Marx, K., Sommer, A., Walter, S. F., Bock, H. G., & Körkel, S. (2018). State and
state of charge estimation for a latent heat storage. Control Engineering Practice, 72, 151-166. doi: 10.1016/j.conengprac.2017.11.006.
The content is identical to the published paper but without the final typesetting by the publisher.
1. Introduction
Thermal energy storages (TES) have a great potential in energy conserva-

tion in the building and industrial sector [1]. Their integration in thermal and
electric/ thermal systems can optimize the use of renewable or industrial waste
energy by peak shaving and shifting strategies leading to a more rational use of
energy and reduction of CO2 emissions [2].
TES using phase change materials (PCM) have been extensively studied re-
cently. Such latent heat TES (LHTES) combine a high energy density with
the advantage of the isothermal nature of the storage process. Applications
of LHTES in buildings primarily aim at improving the performance of space
heating and/or cooling and domestic hot water generation. Typically LHTES
are integrated in sorption systems and seasonal storages, solar collectors, wa-
ter tanks, packed beds, and duct networks [3, 4, 1]. However, PCM might
also be directly integrated into evaporators or condensers, e.g. in domestic re-
frigerators [5]. Industrial applications integrate LHTES in concentrated solar
thermal power plants, and air-conditioning, cold production, and refrigeration
units [2, 6, 7]. Another industrial use case is the exploitation of industrial (in-
termittent) waste heat by transportation using mobile latent heat accumulators
[2].
Under the different methods for storing the thermal latent heat, the change
from solid to liquid has been most widely studied and used [1, 6]. Typical
solid/liquid PCM suitable for low temperatures are e.g. paraffins, fatty acids
and salt hydrates [1]. For high temperature (> 100 ◦ C) applications usage of
inorganic salts, salt eutectic compounds, metal alloys and metallic eutectics is
reported [6]. Each PCM presents its own advantages (e.g. high specific latent
heat and thermal stability, low volumetric expansion during phase change, low
costs) and limitations (such as subcooling and hysteresis effects during phase
change), so its selection has to be done based on the requirements of the respec-
tive application [8, 6, 1]. However, the most serious drawback of solid/liquid
PCM certainly is the poor thermal conductivity. Various techniques for improv-
ing the heat transfer (e.g. encapsulation of PCM or fins and extended surfaces)
have been proposed [8, 6].
Control of thermal systems with storages: Because systems that use energy stor-
ages and the storage itself are inherently transient, effective operating strategies
for dynamic heat integration are needed, such as optimal thermal regulation
functions under transient operating conditions [9, 10]. Model based techniques
and advanced control (such as model predictive control) of TES can significantly
improve the management of dynamic power demands and intensify intermittent
energy services, e.g. in solar thermal systems, district energy systems and build-
ing applications [9, 11]. There is also a growing interest in using TES systems
for demand-side management in the building sector [3]. The idea is to couple
TES with electrically driven heating and cooling systems and to optimally man-
age electrical loads by means of control, adaptation and/or enhancement of the
comparatively large heating, cooling and hot water demands [3].
2
Limitations of current models: Dynamical models of TES and related equip-
ment that can be packaged into commercial or open source software, as well
as accurate models for model-based observers and control can be considered
a key technology for plant and network design, demand side management and
advanced system control [9]. During the last decade, much effort was made
to integrate PCM models into commercial software packages for domestic and
building applications [12]. It seems that corresponding LHTES models are de-
veloped mainly in view of an optimal storage design and integration. They are
frequently based on elaborate computational fluid dynamics and multi-physics
tools, which leads to extremely time-consuming simulations [12]. Still, develop-
ments for model-based observers and controllers seem rather scarce.
State of charge (SOC) of latent thermal storages: Undoubtedly, monitoring of

SOC of a LHTES is a key issue for optimal management of TES and TES sys-
tems. For LHTES with solid/liquid PCM, it seems most reasonable to define the
SOC as the mean liquid phase fraction (or mean melting fraction) of the storage.
The SOC then quantifies the extent to which a latent storage is charged rela-
tive to storable latent heat. This information is crucial for estimating how long
the storage can continue to supply or store latent energy at a given heat flow
and/or temperature level. Moreover, estimation of storage-internal temperature
fields might allow more accurate and tailored charging and discharging opera-
tion strategies, e.g. for realizing complex cyclic and partial loading/unloading
scenarios.
A thorough literature review on the use of the term ‘SOC’ of a latent storage
revealed that this term is rarely used and not unambiguously defined. However,
the terms ‘solidified mass fraction’ or ‘liquid mass fraction’ are frequently used,
describing essentially the same as SOC. Oró et al. [13] use solidified mass fraction
besides other parameters (such as time for complete solidification and number
of transfer units) to characterize the state of an TES system and numerically
investigate the effects of storage-internal design, material and operational pa-
rameters. Similarly, Tiari et al. [14] use ‘liquid fractions’ in a recent numerical
study for a finned heat pipe-assisted TES system. Strumpf et al. [15] use SOC
to denote the melted fraction of PCM for discussing characteristics of a solar re-
ceiver which incorporates a thermal storage with PCM, to be used for electrical
power production for Space Station. SOC is furthermore used in characterizing
the thin layer of PCM for the overheating protection of facade integrated solar
thermal collectors [16] and for characterization of shell and tube (S&T) latent
storage [17]. Vaca Jiménez [18] uses two different definitions for the SOC of a
refrigerator for charge and for discharge (both considering also sensible heat).
They are defined as the ratios of available to maximal energies calculated from
actual and reference temperatures, respectively. Cuneo et al. [19] use the term
‘SOC’ to denote the calculated sensible enthalpy content in a (sensible) stratified
water thermal storage tank, obtained from discrete temperature measurements.
For solar heat receivers employing encapsulated PCM, Hall et al. [20] define the
‘thermal SOC’ as ratio of instantaneously available power to minimum required
power. Henze et al. [21] denote by SOC the ice level in the tank of a chilled-
3
water or ice storage in a cooling system and use this quantity in supervisory
model predictive control for optimal storage charging and discharging.
Measurement based approaches for estimating SOC: Reports on measurement

based approaches for monitoring the phase transition or the phase fractions in
a closed storage during operation are scarce. A possible approach is tracing the
temporal solid/liquid interface through direct visual measurements, see e.g. [22]
and references therein. Henze et al. [21] use direct measurement of the ice level
in a storage tank. Steinmaurer et al. [23] discuss approaches for determining
the SOC and enthalpy of a PCM in closed storage systems. They propose to
apply pressure sensors to measure volumetric changes during phase transition,
which can be feasible for materials with high density differences between liquid
and solid phase (e.g. paraffin). The authors also propose local techniques which
use mechanically induced acoustic waves that propagate through the PCM and
measure the change in propagation speed and the dampening of the signal [23].
Similarly, Bauer et al. [24] propose local electric resistance measurement tech-
niques. Steinmaurer et al. [23] argue that local temperature measurements are
generally not feasible to gather the SOC as this implies high uncertainties and
requires a large number of sensors to get a spatial resolution. Lastly, measure-
ment of heat flux through the storage heat exchanger is also discussed [23].
This approach requires integration of charging and discharging energies while
considering thermal losses. It shows severe disadvantages, since errors in the
estimated thermal losses accumulate with increased application time [23].
State observer applied to thermal storages: Up to the best of the authors knowl-
edge, state and/or parameter observers have not been applied to LHTES until
now. This is possibly due to the complexity of corresponding LHTES models,
being distributed parameter models with (storage-internal) temperature and
phase fractions fields. Kreuzinger et al. [25] present the design of two differ-
ent state observers for a (sensible) stratified water tank for domestic hot water
storage. The observer is used for reconstruction of storage internal time-varying
vertical temperature profiles from a few measurements. The tank is described
by an energy balance equation model which consists of a 1D quasi-linear partial
differential equation (PDE). Kreuzinger et al. [25] adopt a so called late lump-
ing and an early lumping observer design approach. In the former a spatially
weighted correction function is injected in the energy balance equation and the
correction gain is designed based on a physical interpretation. The corrected
model is then numerically solved. In the latter the energy balance equation is
first transformed into an ordinary differential equation (ODE) system by spatial
discretization and an Unscented Kalman filter is used as observer. It is found
that both observer give comparable results in terms of convergence and accuracy.
Considering an optimal placement of storage internal temperature sensors, it is
found that at least three temperature sensors are necessary to attain satisfying
estimation results during all operational modes.
4
This contribution: This paper proposes a definition of the SOC of LHTES with
solid/liquid PCM. A reduced model and a state observer are developed for the
dynamic reconstruction of 2D temperature fields and estimation of SOC of a
lab-scale S&T LHTES. Theoretical and experimental results from testing the
observer are presented.
The paper is organized as follows: Section 2 briefly presents the storage
container design and instrumentation. The energy balance equations for mod-
eling the LHTES are formulated according to previous works. This mathe-
matical model consists of three coupled PDEs and is referred to as detailed
model (dMod). dMod serves in this work as a reference and is used for in
silico studies. In addition, a (physical) simplification of dMod is presented,
which yields a reduced model (rMod) consisting of 12 state variables originat-
ing from the discretization of two coupled PDEs. It is developed with the aim
to design a nonlinear state observer, following the early lumping approach as
described by Kreuzinger et al. [25]. In addition, the definition for the SOC of
LHTES with solid/liquid PCM is given and and its computation is discussed.
Section 3 discusses transformation of rMod into a low order nonlinear ODE
by efficient (coarse) spatial discretization. A collocation scheme is derived for
the 2D heat conduction problem in the PCM cylindrical shell which provides
convenient quadrature formulas to derive the integral storage property SOC.
Section 4 gives details on the numerical implementation. In section 5, predic-
tions from simulations with rMod and dMod are compared to investigate the
influence of the model reduction on the accuracy of predictions. In section 6,
the design and implementation of a continuous-discrete extended Kalman filter
(EKF) is presented and notes on the observability of rMod are given. Accu-
racy and convergence of the EKF are studied using data from simulations with
rMod and experimental data. In section 7, in silico studies are performed to
compare the accuracy of SOC estimates from 1D temperature measurements
only and from 2D temperature information generated by the EKF. Moreover,
experimental results for SOC estimates using the EKF are discussed. Section 8
gives conclusions and directions for future work.
2. Shell and tube 2D LHTES Model
2.1. Container/storage configuration

PCM are typically placed in rectangular or cylindrical containers [26], the
latter minimizing heat losses [27]. The most frequently used and analysed config-
uration is the S&T design with the PCM on the shell side and the heat transfer
fluid (HTF) circulating inside multiple tubes [27, 26], as it offers comparatively
low manufacturing complexity [6] and has shown a superior performance com-
pared to a single tube design [27]. Besides the selection of a suitable PCM,
the cost effectiveness of a solid/liquid LHTES depends mainly on appropriate
measures to increase the heat transfer between HTF and PCM, as the system
performance is limited by the poor thermal conductivity of the PCM [8]. The
following heat transfer enhancement techniques have been employed: extending
5
the heat transfer surface using PCM capsules (normally for low temperature
materials such as paraffins), using PCM embedded porous matrices (metallic
matrix or a matrix made of a naturally porous material such as graphite), dis-
persion of highly conductive particles within the PCM, employing multiple PCM
or using fins [8, 6], employing intermediate heat transfer media or heat pipes
[22].
In this contribution, a lab-scale S&T heat storage with a (latent) storage
capacity of about 10 - 15 kWh is considered. This storage has been described
previously in Zauner et al. [28], Barz et al. [17]. The storage uses a high density
polyethylene (Rigidex HD6070EA) as PCM with phase change temperatures
between 120 and 135 ◦ C. It is a rectangular container with 72 finned tubes
arranged in parallel as a tube bundle. Each tube has 2.7 m length and an
inner and outer diameter of 13.5 and 16.5 mm, respectively. The mean distance
between tubes is 42.8 mm (center to center). The HTF flows through the tubes
and the PCM is located in a shell around the tubes, see figure 1. The container
is equipped with four internal PCM temperature sensors (one-dimensionally
arranged in axial direction), two external HTF temperature sensors, and a HTF
mass flow sensor.
Figure 1: Shell and tube storage design with the HTF flowing through the tubes and the
PCM at the shell side, for detailed information see Zauner et al. [28], Barz et al. [17]. The
figure shows a section of the tube bundle in the rectangular container. ‘PCM 1’ to ‘PCM 4’
indicate positions of four internal PCM temperature sensors. These four sensors are arranged
in one line in axial direction exactly in the middle between neighboring tubes and fins, at 0.1,
0.87, 1.64, 2.4 m. ‘HTF in’ and ‘HTF out’ indicate the position of external HTF temperature
sensors. ‘HTF flow’ is the total HTF mass flow through all tubes.
2.2. Energy balances for the detailed single tube model (dMod)
The following equations are referred to as the detailed model (dMod), serving
as a reference in this work. It is used for in silico studies (when measurements
are not available) representing the real system. In addition, physical simplifica-
tions of dMod yield the model equations of a reduced model (rMod). rMod is
developed with the aim to design a nonlinear state observer.
The tube bundle consists of 72 tubes. For each tube the same inflow and out-
flow conditions are assumed. In addition, it is assumed that all tubes are equal
6
and no heat is exchanged between adjacent tubes. This common assumption
is justified in the present setting due to the (outer) thermal insulation and the
relative abundance of 72 tubes. Based on these assumptions, only one tube with
a surrounding PCM shell needs to be considered for the modeling (single tube
model). To this end, an average fluid flow and effective boundary conditions
are used. The hexagonal shape of the PCM shell is approximated by a cylinder
(from rout to rend ), see figure 2. Finally, a zero heat flux over the outer PCM
shell radius (r = rend ) is enforced. Note that rend marks the radial position of
the installed storage internal PCM temperature sensors (PCM 1-4).
As proposed in previous works [29], energy balances are formulated for three
domains: the HTF, the tube wall, and the PCM (figure 3). The monitoring
of temperature gradients in S&T and pipe TES indicate an essentially two-
dimensional heat transfer in both systems [27]. As largest temperature gradients
are to be expected in the PCM, this domain is modelled in 2D axial and radial
directions. The tube wall is modeled in 1D axial direction with isothermal
conditions in radial direction, axial heat conduction and heat transfer to and
from the PCM and HTF domain. The HTF is modeled in 1D axial direction
with a forced convective flow neglecting heat conduction in axial and radial
direction. It is assumed that the HTF is uniformly distributed to the tube.
The outer storage wall is endowed with a thermal insulation; heat losses to the
surroundings are neglected.
cut A-A
A PCM 1 to PCM 1 to PCM 4
PCM 4
A
L
d
en
t
r
rou
HTF rin
tube
PCM
Figure 2: Shell and tube storage design with the HTF flowing through the tubes and the
PCM at the shell side. Top: Section of the tube bundle and position of the four internal PCM
temperature sensors (‘PCM 1’, to ‘PCM 4’). These four sensors are arranged in one line in
axial direction exactly in the middle between neighboring the tubes. Bottom: Single tube and
modeling domains.
7
PCM
Tube
HTF
Figure 3: Axial and radial modeling domains and differential volume elements of the single
tube model. See table 2 for a list of all variables and their notations.
The mathematical modeling of the conductive heat transfer in the PCM

employs the so-called weak formulation approach, in which the governing heat
equation reads as a single phase equation without explicitly treating the mov-
ing solid/liquid interface. A mushy transition zone between the two phases is
considered, and the apparent specific heat capacity c̃P (T ) is used to model heat
transfer with phase change [30, 31, 32, 33]. Using apparent specific heat capac-
ity (or alternatively, apparent specific enthalpy) is especially useful for technical
grade materials and mixtures that exhibit phase change temperature spans of 10
K or more [33]. This approach is based on the assumption that the PCM mate-
rial undergoing a phase change can be characterized by local (solid/liquid) phase
fractions. In other words, there is no exact temperature (nor an exact location)
but rather a temperature range (or a mushy zone) where the phase transition
occurs. Thus, ’mixed’ thermo-physical properties are defined as functions of the
phase fraction. Consequently, the heat equation can be applied over the entire
region of the domain where the phase change takes place [34, 31, 30]. This pop-
ular approach allows straightforward numeric implementations applying fixed
grid discretization of the PCM domain.
The model equations and boundary conditions are summarized below. The
set of eqs. (1) to (4) are here referred to as detailed model (dMod). Specific
details can be found in Barz et al. [17]. For a list of variables and their notations
see table 2. Unless otherwise stated, volumetric mass densities, specific heat
capacities and thermal conductivities are modeled as functions of temperature.
Balance equation for heat transfer fluid (subscript H)
∂TH ∂TH 2
ρH cp,H = −vH ρH cp,H − q̇H on 0 ≤ x ≤ L (1a)
∂t ∂x rin
0 in
TH (t, x) t=0 = TH (x), TH (t, x) x=0 = TH (t) (1b)
8
Balance equation for tube wall (subscript W)

∂TW ∂ ∂TW 2(rin q̇H − rout q̇P )
ρW cp,W = λW + 2 − r2 on 0 ≤ x ≤ L (2a)
∂t ∂x ∂x rout in
0 ∂TW (t, x) ∂TW (t, x)

TW (t, x) t=0
= TW (x), = 0, =0 (2b)
∂x x=0 ∂x x=L
Balance equation for phase change material (subscript P)

∂TP 1 ∂ ∂TP ∂ ∂TP rout ≤ r ≤ rend
ρP c̃P = r λP + λP on (3a)
∂t r ∂r ∂r ∂x ∂x 0≤x≤L
TP (t, r, x) t=0
= TP0 (r, x), TP (t, r, x) r=rout
= TW (t, x) (3b)
∂TP (t, r, x) ∂TP (t, r, x) ∂TP (t, r, x)

= 0, = 0, =0
∂r r=rend ∂x x=0 ∂x x=L
The heat flux between HTF and tube wall (q̇H ), as well as tube wall and PCM
(q̇P ) is
q̇H (t, x) = α(t, x) (TH (t, x) − TW (t, x)) (4a)

∂TP (t, r, x)
q̇P (t, x) = −λP (t, r, x) r=rout
· (4b)
∂r r=rout
with α being the heat transfer coefficient function.
2.3. Energy balances for the reduced single tube model (rMod)
rMod is a (physical) simplification of dMod. It contains equations for the
HTF, eq. (1), and the PCM, eq. (3). The balance equation for the tube wall
(eq. (2)) is neglected. Instead, a direct heat exchange between HTF and PCM
0
(q̇H ) is assumed
0

q̇H (t, x) = α(t, x) TH (t, x) − TP (t, r, x)|r=rout (5)
0
This means that q̇H (t, x) in eq. (1a) is replaced by q̇H (t, x). In addition, the
second Dirichlet boundary condition in eq. (3b) is replaced by an alternative
0
condition derived from q̇H (t, x) = q̇P (t, x):
λP (t, r, x) r=rout ∂TP (t, r, x)

TP (t, r, x)|r=rout = · + TH (t, x) (6)
α(t, x) ∂r r=rout
9
2.4. Heat capacity and state of charge
A generic two phase model for the apparent specific heat capacity c̃P of the
PCM is given by a linear superposition of terms for liquid (clp,P ) and solid (csp,P )
heat capacity as well as the latent heat (∆H f ) released or absorbed in the phase
transition region [35]:
∂ξ
c̃P := ξclp,P + (1 − ξ)csp,P + ∆H f (7)
| {z } |∂TP{z }
sensible heat
latent heat
where ξ = ml /(ms + ml ) is the liquid weight (mass) fraction, or phase fraction,

and ms and ml represent relative weight masses of solid and liquid state, respec-
tively. For technical grade polymers as considered in this case study the phase
change does not occur at an exact temperature, but rather within a specific
temperature range. This behavior can be described e.g. by defining dξ/dTP in
eq. (7) as a Gaussian pulse [36]. However, because of the higher shape flexibility
in this contribution an adapted Weibull density function is used (see Figure 4
(top)):
( γ−1 γ
∂ξ(TP ) γ µ−TP µ−TP
:= a a exp − a , TP ≤ µ (8)
∂TP 0, TP > µ
where µ, γ, a are location and shape parameters. Details on the choice of the
probability density function representing the phase transition and determining
appropriate parameter estimates from calorimetric measurement data are given
in Barz et al. [17].
Definition (State of charge of a latent thermal storage with solid/liquid PCM).
The term state of charge (SOC) and the symbol Ξ are used to indicate the extent
to which a LHTES is charged relative to storable latent heat. The SOC Ξ is
calculated as the geometric mean of local phase fraction fields ξ(r, x), where r, x
represent spatial coordinates of the distributed parameter system.
With ξ being a function of TP , see eq. (8), local ξ(r, x) are then functions of
local TP (r, x) and are obtained by integration.
Z TP (r,x)
∂ξ(T )
ξ(r, x) = dT (9)
−∞ ∂Tp
For the Weibull density function in eq. (8) a closed form of the cumulative
density function in eq. (9) exists and this integral can be evaluated analytically
[37]. Figure 4 illustrates this relation for the Weibull density. The state of
charge Ξ is calculated as geometric mean (for the cylinder shell geometry) as:
R L R rend
0 rout
ξ(r, x)r dr dx
Ξ= R L R rend (10)
0 rout
r dr dx
10
The integration of eq. (10) needs to be carried out numerically, see below.
0.15
d9(T )=dTP
Weibull
0.10
0.05
0.00
1.00
0.75
9(T )
0.50
0.25
0.00
100 110 120 130 140 150
T [/ C]
Figure 4: Top: probability density function taken from Barz et al. [17].
Bottom: corresponding cumulative density function describing the phase fraction ξ.
3. Spatial discretization of the reduced single tube model (rMod)
Heat transfer in the balance equation of the PCM, eq. (3), is purely con-
ductive. Collocation techniques have proven to be well suited for the numerical
solution of problems of this type [38, 39]. In contrast, axial heat transfer in
the HTF, eq. (1), is purely convective. Thus, discretization by an upwind finite
differences scheme is proposed which is a well established technique in model-
ing fluid flows [40]. In the following, the application of the collocation method
[41] to the balance equations of the PCM (axial and radial direction), eq. (3),
is discussed. For the treatment of the HTF equations the reader is referred to
Quoilin et al. [40], especially for details on the consideration of flow reversal and
two phase flows.
Collocation is a method of weighted residuals (MWR). For conduction or dif-
fusion dominated problems, collocation gives more accurate results and requires
less computation time compared to finite differences computations [41]. For a
discussion of the accuracy of computed radial temperature profiles in the PCM
shell see Appendix A. In the collocation method, the solution is not derived in
terms of the coefficients in the trial function but in terms of the value of the solu-
tion at the collocation points. The independent spatial variables are discretized
and the approximations for the spatial derivatives are defined based on state
values at discretization points. In doing so, differential equations are reduced
to a set of matrix equations. Compared to other MWR, e.g. Galerkin Finite
Element method, the collocation method can be directly applied to nonlinear
equations as no values of integrals need to be evaluated numerically.
In the orthogonal collocation (OC) method the collocation points are taken
as the roots of orthogonal polynomials (for certain weighting functions) [41,
page 97]. In contrast to OC on finite element methods (OCFE), OC uses single
high-order polynomials covering the full domain rout ≤ r ≤ rend , 0 ≤ x ≤ L.
11
In high order approximations, the choice of collocation points is not crucial,
however, using appropriate weighting functions calculations can be made both
convenient and accurate. Moreover, the polynomials can also be generalized
to planar, cylindrical, or spherical geometries as well as to satisfy boundary
conditions. Furthermore, accurate quadrature formulas exist that allow the
exact calculation of integral properties with only a low number of terms. Details
are given in section 4.3 and the quadrature formulae are derived in Appendix
D similarly as in Finlayson [41]. This is important for the presented case study,
as the primary information desired from the solution is an integrated property,
that is the SOC Ξ.
3.1. Normalization of the independent space variables

To apply the collocation method, the spatial differentials in the balance
equations (1),(2),(3) and the integral term in eq. (10) are first transformed into
a standard form with normalized space variables, which can then be directly
approximated by polynomials [38, 39]. The independent space variables r, x are
substituted by the normalized independent space variables R, X ∈ [0, 1].
r − rout x
R=1− , X= (11)
rend − rout L
Note that in eq. (11), the normalized radial coordinate R is mirrored, with R
going from outer (R = 0) to the inner (R = 1) cylinder shell, see figure 5. It will
be shown later that, using a collocation method, it is useful to have the zero
flux (Neumann) boundary condition at position R = 0 rather than at R = 1.
Figure 5: Original and normalized axial and radial coordinates for the PCM shell domain.
Before normalization of space variables, the diffusion terms in x in eq. (2a)

and in r, x in eq. (3a) are transformed in order to eliminate the temperature
12
dependency of λP = λP (TP (r, x)) from the inner differential by applying the
chain and product rule. For the PCM in eq. (3a) this reads:
2
1 ∂ ∂TP dλP ∂TP 1 ∂ ∂TP
r λP = + λP r (12)
r ∂r ∂r dTP ∂r r ∂r ∂r
| {z } | {z }
=:Φ(r,x) =:Ψ(r,x)
2
∂ ∂TP dλP ∂TP ∂ 2 TP
λP = + λP
∂x ∂x dTP ∂x ∂x2
| {z }
=:Ω(r,x)
3.2. Collocation scheme for the PCM

Polynomials suitable for the collocation method are defined on the normal-
ized coordinates R and X, see figure 5. The boundary conditions for the PCM
model require that the first derivative of TP is zero at R = 0, X = 0 and X = 1,
see eq. (3b). For R = 0 and X = 0, these conditions can be directly considered
using symmetric trial functions in R and X, e.g. polynomials of only even pow-
ers, as proposed by Finlayson [38]. The following polynomials that are functions
of R2 and X 2 are used:
NX
R +2 NX
X +1
TP (R ) =2
di R 2i−2
; 2
TP (X ) = d¯k X 2k−2 (13)
i=1 k=1
NR = 1 and NX = 4 define the number of interior collocation points as well

as the number of trial functions to be added and their degree. The placement
of collocation points is shown in figure 6 as part of the complete discretization
scheme for PCM, tube and HTF.
Figure 6: Discretization scheme for PCM and HTF applied to rMod. Pairs of collocation
points (R1 , X1 ), (R1 , X2 ), (R1 , X3 ), (R1 , X4 ) match the positions of the installed storage-
internal PCM temperature sensors. Balance equations are solved and state variables are
computed at collocation points marked by a cross. Elements of the upwind finite differences
scheme for the HTF have been placed such that HTF state variables are computed at the
same axial positions X1 , X2 , X3 , X4 .
13
In the PCM shell, one interior point in R (0 < R2 < 1) and four interior
points in X (0 < X1 < X2 < X3 < X4 < 1) are chosen. Orthogonal collocation
in radial direction leads to R2 ≈ 0.57735 for the interior collocation point (root
of the third Legendre polynomial, weight function w(R) = 1 − R2 for cylindrical
geometry). Note that in Finlayson [41] the weighting function W = 1 − R2 is
proposed for low order approximations which then determines the exact location
of the interior collocation point R2 = 0.57735 from the orthogonality condition
of the trial function TP (R2 ). An additional collocation point is introduced
at R1 = 0. Non-orthogonal collocation is used in axial direction with four
interior points chosen as the (transformed) position of the sensors, such that the
points (R1 , X1 ), (R1 , X2 ), (R1 , X3 ), (R1 , X4 ) exactly match the position of the
four installed storage-internal PCM temperature sensors. Thus, corresponding
computed state variables can be directly compared to measurements, e.g. for
model validation.
This is crucial because the polynomials in eq. (13) deliver a good approxi-
mation of the solution of the energy balance equations only at the collocation
points. For details see Appendix A. Using only one interior collocation point in
radial direction (quadratic polynomial approximation of the radial temperature
profile), an extrapolation to R = 0 (where temperature sensors are installed)
gives poor results. In contrast, if an additional collocation point is placed ex-
actly at R = 0 (quartic polynomial approximation), the polynomial is forced
to fulfill both, the zero flux (Neumann) boundary condition in eq. (3b) and
the energy balance equation in eq. (3a). This additional point/equation then
greatly improves the prediction of the temperature profile. This issue has also
been discussed in Segall et al. [42].
4. Details on the numerical implementation (rMod)
4.1. Normalized equations

Using eq. (11), the following transformed balance equation and initial and
boundary conditions are obtained for the HTF in eq. (1)
∂TH 1 ∂TH 2 0
ρH cp,H = −vH ρH cp,H − q̇ on 0 ≤ X ≤ 1 (14a)
∂t L ∂X rin H
0 in
TH (t, X) t=0 = TH (X) , TH (t, X) X=0 = TH (t) (14b)
14
For the PCM in eq. (3) one gets (see Appendix B), using transformation in
eq. (12) and δ = rend − rout ,
" 2 #
∂TP 1 dλP ∂TP
ρP c̃P = 2 (15a)
∂t δ dTP ∂R
| {z }
=:Φ(R,X)

λP 1 ∂ ∂TP (rout + δ) 1 ∂TP
+ 2 R −
δ R ∂R ∂R (rout + δ(1 − R)) R ∂R
| {z }
=:Ψ(R,X)
" 2 #
1 dλP ∂TP λP ∂ 2 TP 0≤R≤1
+ 2 + 2 on
L dTP ∂X L ∂X 2 0 ≤ X≤ 1
| {z }
=:Ω(R,X)
TP (t, R, X) t=0
= TP0 (R, X) (15b)
∂TP (t, R, X) ∂TP (t, R, X) ∂TP (t, R, X)
= 0, = 0, =0
∂R R=0 ∂X X=0 ∂X X=1
0
and for the heat flux, q̇H in eq. (5), one gets
0

q̇H (t, X) = α(t, X) TH (t, X) − TP (t, R, X)|R=1 (16)
The boundary condition for rMod in eq. (6) reads
λP (t, R, X) R=1 ∂TP (t, R, X)

TP (t, R, X)|R=1 = − · + TH (t, X) (17)
α(t, X)(rend − rout ) ∂R R=1
The normalization of the integral term in the numerator in eq. (10) gives,
R1
L 0 Γ(X) dX
Ξ= 2 2 )/2 (18a)
L(rend − rout
Z Z
2
1 2
1
Γ(X) = δ + δrout ξ(R, X) dR − δ ξ(R, X)R dR (18b)
0 0
Z TP (R,X)
ξ(T )
ξ(R, X) = dT (18c)
−∞ dT
4.2. Collocation applied to the PCM shell

With NR = 1 and NX = 4, the collocation scheme of the PCM shell gives
eight discrete equations, see figure 6:
∂TP
ρP c̃P = Φ(Rj , Xl ) + Ψ(Rj , Xl ) + Ω(Rj , Xl ) (19)
∂t Rj ,Xl
j = 1, NR + 1 ; l = 1, · · · , NX
15
Note that subscripts j and l are used to denote individual collocation points in
R and X, respectively. At these points the spatial derivatives in Φ(Rj , Xl ) and
Ψ(Rj , Xl ) in eq. (19) (see eq. (15a) for their definition) can be replaced by the
collocation derivative formulas:
!2 NX
!2
R +2
∂TP (2)
= Aji TP (Ri , Xl ) (20a)
∂R Rj ,Xl i=1
NXR +2
1 ∂ ∂TP (2)
R = Bji TP (Ri , Xl ) (20b)
R ∂R ∂R Rj ,Xl i=1
NX
R +2
1 ∂TP (2)
= Ãji TP (Ri , Xl ) (20c)
R ∂R Rj ,Xl i=1
The coefficient matrices A(2) , B(2) ∈ R(NR +2)×(NR +2) are given for cylindrical
geometry and are derived from the collocation points Rj with j = 1, · · · , NR + 2
as defined in Finlayson [38]. The calculation of the coefficients of the adapted
matrix Ã(1) is given in Appendix C. The term (rout + δ)/(rout + δ(1 − R)) in
Ψ(Rj , Xl ) in eq. (15a) is evaluated for Rj .
In the same way, the partial derivatives in Ω(Rj , Xl ) in eq. (19) can be
replaced by
!2 NX
!2
X +1
∂TP (1)
= Alk TP (Rj , Xk ) (21a)
∂X Rj ,Xl
k=1
NX
X +1
∂ 2 TP (1)
= Blk TP (Rj , Xk ) (21b)
∂X 2 Rj ,Xl k=1
The coefficient matrices A(1) , B(1) ∈ R(NX +1)×(NX +1) are given for planar geom-
etry and are derived from the collocation points Xl with l = 1, · · · , NX + 1 as
defined in Finlayson [38].
The following boundary conditions are considered in the discrete form of
the model equations for PCM. The zero flux (Neumann) boundary condition
∂TP (t, R, X)/∂X|X=1 = 0 in eq. (15b) is incorporated using the collocation
formula [42]:
PNX (1)
A TP (Rj , Xk )
TP (Rj , Xm ) = − k=1 mk(1) ; m = NX + 1 (22)
Amm
Equation (22) is an explicit expression which is used to eliminate temperatures
TP (Rj , Xk ) with j = 1, 2 and k = NX + 1 = 5 from eq. (21).
The normalized boundary condition for rMod in eq. (17) is discretized re-
placing the spatial derivative with:
NX
R +1
∂TP (2) (2)
= A(NR +2)i TP (Ri , Xl ) + A(NR +2)(NR +2) TP (RNR +2 , Xl )
∂R RNR +2 ,Xl i=1
(23)
16
Using eq. (23) and with NR + 2 = 3, an explicit expression can be derived which
is used to eliminate temperatures TP (R3 , Xl ) with l = 1, · · · , 4 from eq. (20):
λP (R3 ,Xl ) P2 (2)
α(Xl )TH (Xl ) + rend −rout · i=1 A3,i TP (Ri , Xl )
TP (R3 , Xl ) = (2)
(24)
α(Xl ) − λrend
P (R3 ,Xl )
−rout A3,3
4.3. Evaluation of integrals

Local phase fractions in eq. (18c) are computed from the cumulative Weibull
density function, see NIST/SEMATECH [37]. The integrals in eq. (18b) can be
replaced by
Z 1 NX
R +2
(2)
ξ(R, X)R dR = wj ξ(Rj , X) (25)
0 j=1
Z 1 NX
R +2
(2)
ξ(R, X) dR = w̃j ξ(Rj , X)
0 j=1
where the vector w(2) and the adapted vector w̃(2) are given for cylindrical
geometry, see Appendix D. Finally, the integral over X in eq. (18a) is evaluated
using the vector w(1) for planar geometries.
5. Comparison of detailed and reduced single tube models dMod and

rMod
dMod represents a well accepted model complexity for numerical studies of

the heat transfer in S&T configurations with technical grade solid/liquid phase
change materials [29]. rMod is an approximation of dMod with the aim to
design a nonlinear state observer. First, the energy balance of the tube wall is
neglected (see eqs. (5) and (6)), second, the model order is further reduced by
efficient (coarse) spatial discretization of the HTF and PCM energy balances,
see section 3. In this section 5 it is tested if rMod preserves the original dynamic
behaviour of dMod as much as possible. For doing so, model predictions of rMod
and dMod are compared.
dMod has been fitted to experimental data from thermo-physical measure-
ments and from lab-scale LHTES operation. It has proved to agree well with
experimental results, provided that the phase transition behavior of the PCM
can be described by a unique relationship between temperature and phase frac-
tion (as defined in eqs. (8) and (9)). This holds, e.g., if the phase transition
is not noticeably affected by subcooling, which might lead to hysteresis effects
[17]. The detailed model (dMod) is solved using a fine discretization with 50
elements in axial direction and 7 elements in radial direction, yielding a (sparse)
system of 450 ODEs (for detailed information on dMod see Barz et al. [17]).
17
Table 1: List of variables used in model based and experimental analysis. Normalized positions
as used in figure 6 and eq. (27): X = [0.04, 0.34, 0.66, 0.96], R = [0.00, 0.58])
notation variable symbol

in figures element in rMod meaning
in
HTF in u1 TH HTF temperature at storage inlet
HTF flow u2 ṁH,total total HTF mass flow through storage
HTF out y5 TH (X4 ) HTF temperature at storage outlet
PCM 1-4 y1 , · · · , y4 TP (R1 , X1 ) temperatures at the outer edge of the PCM shell
..
. at four different axial positions
TP (R1 , X4 ) (positions match the installed PCM temperature sensors)
PCM 5-8 x9 , · · · , x12 TP (R2 , X1 ) temperatures inside the PCM shell
..
. at four different axial positions
TP (R2 , X4 ) (no measurements available)
rMod yields a system of 12 continuous-time nonlinear state equations and 5

discrete-time measurement equations:
ẋ(t) = f (x(t), u(t), p) (26)

y(tk ) = h(x(tk ))
For convenience, we omit the explicit notation of the time argument t in the
following. In eq. (26), x ∈ R12 are the states, u ∈ R2 are the inputs, y ∈ R5
are the measurements of the system taken at tk and h simply assigns a subsets
of states to y. Using the notation for discrete states as introduced in section 4,
the vector elements of these variables read:
x = [TH (X1 ), TH (X2 ), TH (X3 ), TH (X4 ), (27)

TP (R1 , X1 ), TP (R1 , X2 ), TP (R1 , X3 ), TP (R1 , X4 ),
TP (R2 , X1 ), TP (R2 , X2 ), TP (R2 , X3 ), TP (R2 , X4 )]T
in
u = [TH , ṁH,total ]T
y = [TP (R1 , X1 ), TP (R1 , X2 ), TP (R1 , X3 ), TP (R1 , X4 ), TH (X4 )]T
In eq. (27), the first 4 elements in x represent the HTF temperatures at different
axial coordinates, and the remaining eight states describe the PCM tempera-
tures at different axial and radial coordinates, see figure 6. The input variable
u consists of the temperature at the HTF storage inlet (HTF in) and the HTF
2
total mass flow (HTF flow), which is calculated as ṁH,total = nT vH π (rin ) ρH ,
with nT being the number of tubes. These input variables (denoted as ‘input’ in
the rest of the paper) are considered as time dependent parameters in the model.
The elements in the measurement vector y are used to compare model predic-
tions to real measurement data. Five corresponding temperature measurements
(denoted as ‘meas’ in the rest of the paper) are available: four storage-internal
PCM temperatures (PCM 1-4) and the HTF temperature at the storage outlet
(HTF out), see figure 1 for the position of sensors and figure 6 for the location
of corresponding discrete states. An overview is given in table 1.
18
In the following a comparison of predictions generated by simulations with
rMod and dMod is given. Figure 7 shows results for the HTF temperature at
the storage outlet y5 (HTF out, first subfigure), and results for four different
storage-internal PCM temperatures y1 , · · · , y4 (PCM 1-4, second subfigure).
The dynamic experiment design (input) is defined by variations in u1 , this is
the HTF storage inlet temperature (HTF in, first subfigure), and in u2 , this is
the HTF total mass flow (HTF flow, third subfigure).
160 HTF in (input)

HTF out (dMod)
HTF out (rMod)
T [/ C]
140
HTF in
HTF out
120
160 PCM 1-4 (dMod)

PCM 1-4 (rMod)
T [/ C]
140
120 PCM 1, 2,3,4
2
_ [kg/sec]
1 HTF flow (input)

m
0
0 61.7 100 150 200 236.8 300
time [min]
Figure 7: Predictions from simulations with the detailed model (dMod, see Barz et al. [17])
and the reduced model (rMod) for a dynamic experiment design (input) defined by variations
in the HTF storage inlet temperature and HTF total mass flow show only minor differences.
It can be seen in figure 7 that the predicted HTF out (dMod) and HTF out
(rMod) (first subfigure), as well as the predicted PCM 1-4 (dMod) and PCM
1-4 (rMod) (second subfigure) agree well. The largest deviations can be found
for fast changes in HTF out (around t = 236.8 min) with a maximal absolute
difference of 3.1 K and an average absolute difference of 0.6 K. For the predicted
PCM 1-4 temperatures the absolute differences have a maximal value of 4.2 K
and an average value of 0.4 K.
Predicted temperature fields in the PCM shell are depicted in figure 8 for
two selected time points (snapshots). Temperature values at collocation points
at r = rend correspond to temperatures PCM 1-4 in figure 7 (the four points
where the storage-internal temperature sensors are installed). Note that the
balance equations eq. (3) are fulfilled only at the collocation points.
19
160 145
155
T [/ C]
T [/ C]
150 140
145 rend =0.021 rend =0.021
140 0.014 135 0.014

2.5 2.5 2
2 1.5 1 rout =0.008 r [m]
1.5 1 0.5 rout =0.008 r [m] 0.5 0
x [m] 0 x [m]
(a) (b)
Figure 8: Predictions from simulations with the reduced model (rMod) for the dynamic ex-
periment design in figure 7. The figure shows temperature fields in the PCM shell at time
61.7 min (a) and time 236.8 min (b).
Temperature values at collocation points (indicated in figure 6) are marked by a star. Con-
tinuous lines in axial and radial direction are obtained from inter-/extrapolation using the
trial functions in eq. (13). Radial temperature differences are relatively small here due to the
operating conditions.
6. The state observer
An overview of currently available observer methods for state estimation of

nonlinear systems can be found e.g. in Rawlings and Bakshi [43]. The idea is to
use feed-back control from available HTF and PCM temperature measurements
(y in eq. (27)) to provide state estimates of not measured storage-internal HTF
and PCM temperatures (xi ∈ / y, with i = 1, · · · , 12 in eq. (27)) and to improve
the predictions generated by simulations with rMod.
rMod has continuous-time state equations and discrete-time measurement
equations, see eq. (26). A continuous-discrete extended Kalman filter [44], also
called hybrid EKF [45], is used to estimate the states from noisy measurements.
To this end the process and the measurement errors in eq. (26) are assumed to
be independent Gaussian white noises with zero mean and covariance matrices
Q and R.
6.1. State observer design and observability

Because of the complexity of the reduced nonlinear model (rMod), observ-
ability is not studied in a mathematical rigorous way. However, convergence of
the filters is studied for a wide operating range of several simulated and real
input signals (HTF inlet temperature and HTF mass flow). Local evaluation of
Kalman rank condition and Hautus test indicate full rank. In addition, conver-
gence is also successfully tested for deviations in the initial conditions between
the observer and the plant, see the following section 6.3 for an illustration.
rMod is a reduced model with 12 states. It has been reduced such that it is
observable for the number of actuators and sensor points in the studied LHTES
(see eq. (27) and table 1). Using additional collocation points, either in axial
or radial direction, or considering the energy balance equation of the tube wall,
20
increases the number of states. It has been found that the resulting system is
then no longer observable, the EKF diverges. The same applies if less sensor
points are considered, i.e. when the number of temperature sensors is reduced.
Validation of the observer can only be performed with measurement data
that has not been fed into the observer. Unfortunately, the lab-scale LHTES
has only been endowed with a minimal number of temperature sensors. As con-
sequence, a validation for xi ∈
/ y, with i = 1, · · · , 12, cannot be done. However,
the accuracy of such estimations is assessed by in silico studies using experi-
mental data generated by simulations with dMod, see Appendix E.
6.2. Numerical implementation

The observer equations are implemented as described in Schneider and Geor-
gakis [45]. The EKF parameters have been selected as follows: diagonal elements
in the 12 × 12 initial error covariance matrix P0 and in the 12 × 12 model noise
covariance matrix Q are 1.11 × 10−2 ; diagonal elements in the 5 × 5 measure-
ment noise matrix R are 1.111 × 10−1 .
All computations have been carried out in Matlab on an Intel R CoreTM 2 Duo
Processor E6600 @2.40-GHz computer with 4 GB RAM. The system Jacobian
is approximated using finite difference approximations. The filter performance
is examined in in silico studies with measurements generated using dMod, and
in experimental studies with measurements from a lab-scale LHTES. All mea-
surements are stored at discrete time steps tk = k∆t with ∆t = 25 sec. In silico
and experimental studies have been both performed using stored data.
∆t are selected such that the resulting time series mimics the important
dynamics of the system states, namely TH and TP . The dynamics of TH depend
on the mass flow, with a rough estimate of the typical time constant lying
between one and five minutes. The dynamics of the PCM temperatures are
affected by the huge changes (around factor 10) in the specific heat capacity
during phase transition. A rough estimate of typical time constants gives 10
to 100 min. The overall computation time of the filter is around 22 sec for a
typical experiment duration of 300 min with k = 720 time steps. Thus, the full
filter computation time is even lower than the sample interval.
6.3. Experimental results

The performance of the state observer is tested for experimental data from
the lab-scale LHTES presented previously [17, 28]. Results are distinguished
using the notations:
• rMod: predictions generated by simulations with rMod,

• EKF-rMod: state estimations generated by the observer.
Figures 9 and 10 compare rMod-based simulations to state estimations gen-
erated by the EKF-rMod observer. As can be seen from the re-initialization
events at 10 and 125 min with poor initial data, the rMod-based simulations
suffer from poor state estimates until a switch from heating to cooling occurs
21
(figure 9(a), top), whereas the EKF-rMod observer rapidly converges to the
actual state (figure 9(b), top). The PCM measurements (figure 9(a) and fig-
ure 9(b), bottom) show similar behaviour.
The same can be observed in figure 10(a). Note that here the simulation
and the observer are both correctly initialized during steady operation (0 to 15
min) with temperature values equal 85 ◦ C. However, figure 10(a) also deserves
special analysis, as after 50 min the HTF flow is stopped. The measured external
temperatures HTF in and HTF out show a decrease (from 50 to 100 min). This
can be explained by heat losses to the surroundings, that have been neglected in
the modeling, especially at the container inlet and outlet tubes and baffles. It
can also be seen, that the measured PCM temperature PCM 1 follows HTF in.
Note that, the sensor position of PCM 1 is close to the container inlet with x=0.1
m, see also figure 1. Interestingly enough, the situation changes after 100 min
and the decrease in measured HTF and PCM temperatures now continues on a
slower pace. This can be explained by the fact, that PCM temperatures have
reached the melting temperature range, (see temperatures and corresponding
phase fractions in figure 4), where the latent heat is released.
It is important to note that heat losses are not considered in rMod. Ac-
cordingly, in figure 10(a), after 50 min, the predictions (rMod) for internal and
external temperatures are far off the actual values. Matters are quite different
for state estimations (EKF-rMod), see figures 9(b), 10(b). These state estima-
tions converge much faster (figure 9(b)) and generally agree very well, even if
the HTF flow is stopped (figure 10(b) after 50 min). It can be concluded that
the state observer is able to robustly track real process data even for complex
operating scenarios and conditions neglected in the model.
Note that differences between measurements (meas) and predictions (rMod),
and, measurements (meas) and state estimations (EKF-rMod) are not quanti-
tatively compared, as the tuning parameters of the EKF-rMod observer may be
chosen to either put more (or less) emphasis on actual measurements.
22
160 HTF in (input)
init HTF out (meas)
140
T [/ C]
HTF out (rMod)

120
HTF in
HTF out
100
reinit
80
160 PCM 1-4 (meas)
init PCM 1-4 (rMod)
140
T [/ C]
120
PCM 1,2,3,4
100
reinit
80
0 50 100 150 200 250 300
time [min]
(a)
160 HTF in (input)

init HTF out (meas)
140
T [/ C]
HTF out (EKF-rMod)

120
HTF in
HTF out
100
reinit
80
160 PCM 1-4 (meas)
init PCM 1-4 (EKF-rMod)
140
T [/ C]
120
PCM 1,2,3,4
100
reinit
80
0 50 100 150 200 250 300
time [min]
(b)
Figure 9: (a) Predictions generated by simulations (rMod), (b) state estimations gener-
ated by the observer (EKF-rMod), and corresponding experimental data (meas) for dy-
namic experiment design (input) of HTF inlet temperature and a constant total mass flow of
ṁH,total = 1.09 kg/sec.
23
HTF in (input)
160 HTF out (meas)
HTF out (rMod)
140
T [/ C]
120 HTF out
HTF in
100
80
160
140
T [/ C]
120
PCM 1,2,3,4 PCM 1-4 (meas)
100
PCM 1-4 (rMod)
80
_ [kg/sec]
1
HTF flow (input)
0.5
m
0
0 20 40 60 80 100 120 140 160
time [min]
(a)
HTF in (input)
160 HTF out (meas)
HTF out (EKF-rMod)
140
T [/ C]
120 HTF out

HTF in
100
80
160
140
T [/ C]
120
PCM 1,2,3,4 PCM 1-4 (meas)
100
PCM 1-4 (EKF-rMod)
80
_ [kg/sec]
1
HTF flow (input)
0.5
m
0
0 20 40 60 80 100 120 140 160
time [min]
(b)
Figure 10: (a) Predictions generated by simulations (rMod), (b) state estimations gener-
ated by the observer (EKF-rMod), and corresponding experimental data (meas) for dynamic
experiment design (input) of HTF inlet temperature and total mass flow. Measurement arte-
facts/noise around time t = 70 and t = 100 min do not interfere with the EKF-rMod observer.
24
7. Prediction of the state of charge (SOC)
The SOC Ξ of the storage is calculated as the geometric mean of the phase
fraction field ξ(r, x) in the cylindrical PCM shell, see eq. (10). Figure 11 shows
an example of a temperature field TP (r, x) and corresponding local phase frac-
tion field ξ(r, x) calculated by eq. (9).
150
1
140
T [/ C]
9 [-]
0.5
130
120 0
0 rout =0.008 0 rout =0.008
0.5 0.5
1 0.014 1 0.014
1.5 1.5
2 2
x [m] 2.5 rend =0.021 x [m] 2.5 rend =0.021
r [m] r [m]
(a) (b)
Figure 11: Example of a local temperature field TP (r, x) in the PCM shell (a) and correspond-
ing local phase fraction field ξ(r, x) (b). Temperature TP and phase fraction ξ at collocation
points (indicated in figure 6) are marked by a star. TP and ξ fields are given for the trial
functions. Note that the values of TP and ξ fulfill the energy balances at the collocation points
only.
7.1. In silico studies

A direct evaluation of the accuracy of Ξ obtained from state estimations
(EKF-rMod) is impossible, because the true SOC cannot be measured. However,
in silico studies have been performed, in which the real process is represented
by dMod, and in which the SOC is calculated from simulations with dMod.
Two approaches are assessed regarding their accuracy of estimating the SOC.
The first approach uses data from available measurements only (direct mapping,
dirMap), the second approach uses information generated by the observer (EKF-
rMod). Accordingly, the following methods are considered to obtain the SOC
Ξ:
• dMod: using complete temperature information from dMod, i.e. calcula-
tion of a phase fraction field ξ(r, x) and its geometric mean Ξ from the
simulated local temperature field TP (r, x) (for in silico studies only),
• dirMap: using temperature information from installed PCM temperature
sensors only, i.e. calculation of geometric mean Ξ of the four phase frac-
tions ξ1 , · · · , ξ4 that are obtained from the corresponding temperature
measurements TP at PCM 1-4,
• EKF-rMod: using temperature information generated by the observer, i.e.
calculation of a phase fraction field ξ(r, x) and its geometric mean Ξ from
state estimations of the local temperature field TP (r, x).
25
100
75
SOC [%]
50
25
1 # PCM shell thickness
0
100
75
SOC [%]
50
25
0
100
dMod
75 dirMap
SOC [%]
EKF-rMod
50
25
0
0 50 100 150 200 250 300
time [min]
(a)
1# PCM shell thickness

140
T [/ C]
TP (rout )!
TP (r) at
135 A
t =177 min
TP (rend )
130
140 dMod
T [/ C]
rMod
135
A
130 meas *
140
T [/ C]
135
130
150 177 200 rout rend 2rend 3rend
time [min] r [m]
(b)
Figure 12: In silico studies for three storage designs with increased PCM shell thickness
and for the dynamic experiment design shown in figure 7. The real process is represented by
simulations with dMod. (a) shows results for SOC, (b, left) shows different PCM temperatures
from dMod, between 150 and 200 min, at the axial position of the second installed temperature
sensor PCM 2; (b, right) shows corresponding snapshot radial temperature profiles at 177 min.
The outer edge of the PCM shell is defined by rend , 2rend , 3rend , respectively.
SOC is calculated from available local PCM temperature information: (dMod) uses PCM
temperature fields from the detailed model dMod; (dirMap) uses measurements from simulated
four PCM temperature sensors; (EKF-rMod) uses estimated PCM temperature fields from the
observer.
Figure 12 shows the results for three storage designs with increased PCM
shell thickness. Calculations relying only on measurement data (dirMap) lead
to poor results, see figure 12a. It can be seen that the error increases with
26
increasing shell thickness. For the storage design used in this work (1× PCM
shell thickness) a time delay between calculations from ‘dMod’ and ‘dirMap’
of about 4-5 min is observed. This time delay reaches 15-20 min for 3× PCM
shell thickness. This dramatic increase is to be expected, as the volume of the
PCM increases quadratically with increasing shell thickness due to the cylindri-
cal shell geometry. Accordingly, temperature measurements from the installed
four internal temperature sensors at the outer edge of the PCM shell do not
represent the radial temperature fields and thus the dynamic response times
to temperature changes at the PCM inner shell increase significantly, see also
figure 12b. As a result, absolute SOC values are not reflected correctly and
changes in the SOC are recorded with a relatively large time delay.
In contrast, the observer-based estimation of the SOC Ξ gives correct results
(compare results for ‘EKF-rMod’ and ‘dMod’).
7.2. Experimental studies

Figure 13 shows results for the calculated SOC for real experimental data
from the lab-scale LHTES presented previously [17, 28]. The two different
dynamic experiment designs from figures 9 and 10 are considered.
Similarly as for the results in figure 12, there exists a (relatively small) time
delay between calculated SOC from ‘EKF-rMod’ and from ‘dirMap’ of about
4-5 min, see figure 13 (top). Interestingly enough, in figure 13 (bottom), besides
the time delay, also relatively large differences between the SOC calculated from
‘EKF-rMod’ and from ‘dirMap’ can be observed. Remember that, after 50 min,
the HTF flow is stopped and the storage cools down due to heat losses to
the surroundings, see figure 10(b). The installed internal PCM temperature
sensors detect temperature drops only with a large time delay. As a result, the
calculated SOC (dirMap) hardly changes after 60 min. In contrast, the EKF-
rMod observer detects these heat losses indirectly based on information from
external temperature measurements at the storage inlet (HTF in) and outlet
(HTF out). Internal PCM temperature values are corrected by the observer
and the SOC decreases, see also figure 10(b).
In figure 13 (bottom), it can also be seen that measurement noise directly
affects SOC values generated by dirMod, whilst SOC values generated by EKF-
rMod are filtered and thus not affected.
This fact highlights again the superior performance of the observer based
monitoring of the SOC (EKF-rMod) compared to the approach based on avail-
able information from installed temperature sensors only (dirMap).
27
100 dirMap
EKF-rMod
75
SOC [%]
50
25
0
0 50 100 150 200 250 300
time [min]
100
75
SOC [%]
50 dirMap
HTF flow is stopped EKF-rMod
25
0
0 20 40 60 80 100 120 140 160
time [min]
Figure 13: Computation of SOC of the storage for the dynamic experiment designs shown in
figure 9 (here at top) and in figure 10 (here at bottom).
The SOC is calculated from available local PCM temperature information: (dirMap) uses
measurements from installed four PCM temperature sensors; (EKF-rMod) uses estimated
PCM temperature fields from the observer.
In dynamic charging and discharging (top, cf. figure 9), both the EKF-rMod observer and
the direct mapping approach deliver comparable results. For the load/stop scenario (bottom,
cf. figure 10), the EKF-rMOD observer delivers more convincing results.
8. Conclusions
In this paper, a reduced physical 2D dynamic model of a lab-scale thermal

storage with solid/liquid phase change material is developed. The storage inter-
nals follow the most common S&T design, with the PCM at the shell side and the
HTF flowing through the pipes in the center. The developed model has proven
to be sufficiently accurate and to meet the requirements for efficient system sim-
ulations and for the use with model-based observers and control methods. The
energy balance equation model consists of two coupled PDEs, one for the HTF
in a tube and one for the surrounding PCM shell. Space discretization leads to
a small system of 12 nonlinear ODEs which can be treated by classical solution
strategies, e.g. initial value solvers.
Based on this model, the design and application of a novel state observer for
the reconstruction of transient spatial temperature fields inside the storage is
described. Moreover, to the best of the authors knowledge, the first experimental
application of a state observer to monitor the state of charge (SOC) of a LHTES
is presented.
A collocation scheme is proposed for the numeric solution of the 2D heat
conduction equation of the PCM cylindrical shell. To account for a zero flux
condition at the outer shell diameter, a transformation of the governing PDE
by normalization is derived and adapted coefficient matrices are used in the
collocation scheme. The implementation is not straightforward and, to the best
of the authors knowledge, has not been presented elsewhere. It can be easily
adopted for similar heat conduction and mass diffusion problems in cylinders
28
and cylindrical shells with zero flux and/or symmetry boundary condition at the
outer shell diameter. The applied collocation method yields excellent predictions
already for low order polynomial approximations (in this case study quartic
polynomials are used) of the states and derivatives at collocation points as well
as of integral properties of state profiles (using quadrature formulas). It is shown
how these integral properties can be used to compute geometric mean values,
i.e. mean phase fractions and the SOC of the LHTES.
In this contribution, the discretization grid (collocation points and finite
element boundaries) has been chosen such that corresponding discrete state vari-
ables match the position of installed storage-internal PCM temperature sensors.
Installation of these sensors gives valuable insights and is practiced at least in
lab-scale storages, see e.g. [46]. The optimal number and location of temperature
sensors and the optimal discretization scheme for alternative storage geometries
may be interesting research topics as well as improving the observability by sub-
stituting the EKF by a moving horizon estimator (MHE) that uses more than
only the most recent measurement (see, e.g., the real-time MHE for simultane-
ous estimation of states and parameters in Kühl et al. [47]). Moreover, in silico
studies have been performed to validate observer estimations for non-measured
states (temperatures inside the PCM shell). Thus, further completion could be
an experimental validation by use of additional sensors.
In this contribution, the SOC of an LHTES is defined as the mean of local
phase fraction fields in the solid/liquid PCM in the storage. A simple algebraic
phase transition model is used which is an unambiguous assignment/mapping
of the phase fraction to the PCM temperature. This model can be derived
from PCM thermo-physical material data. However, more complex modeling
is conceivable, e.g. an ambiguous assignment/mapping which could result, e.g.,
from the consideration of subcooling or kinetics in the phase transition. Both
phenomena are of practical relevance, though most often unwanted [7, 6, 46, 48].
Acknowledgements
This work was partly funded by the Austrian Research Funding Association
(FFG) within the programme Bridge in the project ‘modELTES’ (project No.
851262). Support within the Advanced Investigator Grant MOBOCON (grant
agreement No. 291458) of the European Research Council is gratefully acknowl-
edged. The authors want to thank an anonymous reviewer for his thoughtful
and constructive comments on an earlier version of this manuscript.
29
Acronyms and Notation
Table 2: List of variables used for the LHTES model.
symbol extent ∗ units meaning

a — — shape parameter of the Weibull density function
α a, t W K −1 m−2 heat transfer coefficient function
cp,H a, t J K −1 kg −1 specific heat capacity of HTF
cp,W a, t J K −1 kg −1 specific heat capacity of tube wall
cp,P a, r, t J K −1 kg −1 specific heat capacity of PCM
c̃P a, r, t J K −1 kg −1 apparent specific heat capacity of PCM
δ — m PCM shell thickness in radial direction
∆H f — J kg −1 specific latent heat of PCM
γ — — shape parameter of the Weibull density function
L — m storage’s physical length in axial direction
λP a, r, t W m−1 K −1 thermal conductivity of PCM
λW a, t W m−1 K −1 thermal conductivity of tube wall
ṁH t kg s−1 HTF mass flux
µ — — location parameter of the Weibull density function
nT — — HTF tube count in thermal storage
NR , NX — — number of collocation points in radial/axial direction
q̇H a, t W m−2 heat flux density from HTF into tube wall
q̇P a, t W m−2 heat flux density from tube wall into PCM
r, R r m radial coordinate (R: normalized)
rend — m combined PCM shell and HTF tube radius
rin — m inner HTF tube radius
rout — m outer HTF tube radius
ρH — kg m−3 volumetric mass density of HTF
ρW — kg m−3 volumetric mass density of tube wall
ρP — kg m−3 volumetric mass density of PCM
◦
TH a, t C HTF temperature
0 ◦
TH a C HTF initial temperature profile
in ◦
TH t C HTF inlet temperature
◦
TP a, r, t C PCM temperature
◦
TP0 a, r C PCM initial temperature field
◦
TW a, t C tube wall temperature
vH t m s−1 velocity of HTF
x, X a m axial coordinate (X: normalized)
ξ a, r, t — local (solid) phase fraction field
Ξ t — state of charge of thermal storage
∗
extent may be axial, r adial, and temporal
Table 3: List of acronyms
acronym meaning acronym meaning

dirMap direct mapping OCFE OC on finite elements
dMod detailed model ODE ordinary differential equation
EKF extended Kalman filter PCM phase change material
HTF heat transfer fluid PDE partial differential equation
LHTES latent heat thermal energy storage rMod reduced model
MHE moving horizon estimator SOC state of charge
MWR method of weighted residuals S&T shell and tube
OC orthogonal collocation TES thermal energy storage
30
Appendix A. Discretization error
Figure A.1 shows the discretization error for OC schemes differing in the
number and location of collocation points for the one-dimensional nonlinear
dynamic heat transfer problem in radial direction in the PCM cylinder shell.
Numeric results are compared to the benchmark solution computed by a finite
differences scheme with 50 elements, see Barz et al. [17] for details. In fig-
ure A.1(a) results obtained by using a single interior collocation point (which
corresponds to one equation to be solved) are shown. The solutions computed at
the interior collocation point show relatively high errors. Moreover, the extra-
polation to the outer edge of the cylinder shell at rend , where the temperature
sensors are installed, yields poor results. In contrast, figure A.1(b) shows that
the results are much better when using an additional collocation point placed at
rend . The same holds true for the geometric mean temperatures, Tmean , shown
on the right of figure A.1. Their values are obtained from the evaluation of the
integral term in eq. (18b) using the quadrature formulas in eq. (25).
The absolute errors listed in table A.1 confirm the above findings. It can
therefore be concluded that the value of the integral property SOC Ξ can
also be approximated with reasonable accuracy using the OC scheme in fig-
ure A.1 (b). Moreover, the OC scheme generally outperforms the FD scheme in
terms of the accuracy of predictions (compare results with the same number of
points/equations to be solved in table A.1).
Table A.1: Absolute errors in the geometric mean Tmean in K for different discretization
schemes: Orthogonal Collocation (OC) and Finite Differences (FD). Computed errors refer
to deviations from geometric mean values of FD benchmark trajectory (FD with 50 elements,
shown in figure A.1). The total number of interior (int.) and boundary (bound.) points
indicate the number of equations to be solved.
points abs. error [K] at

int. bound. t1 t2 t3
OC 1 - 2.4 2.9 2.8 see figure A.1 (a)
OC 1 1 1.4 0.7 0.2 see figure A.1 (b)
OC 2 1 0.4 0.2 0.0 see figure A.1 (c)
OC 3 1 0.2 0.0 0.0 see figure A.1 (d)
OC 4 1 0.1 0.0 0.0
FD 1 - 5.5 3.4 0.8
FD 2 - 2.7 1.5 0.3
FD 3 - 1.8 0.8 0.2
FD 4 - 1.3 0.4 0.1
FD 5 - 1.0 0.3 0.1
31
140 T as f (r) Tmean 140 T as f (r) Tmean
t3
t3
130 130
t2
T [/ C]
T [/ C]
t2
120 120
t1
t1
OC (collocation point) OC (collocation point)

OC (interpolation) OC (interpolation)
110 110
FD with 50 elements FD with 50 elements
rout = 0.00825 0.013813 rend = 0.021412 rout = 0.00825 0.013813 rend = 0.021412
r [m] r [m]
(a) (b)
140 T as f (r) Tmean 140 T as f (r) Tmean

t3 t3
130 130
T [/ C]
T [/ C]
t2 t2
120 t1 120 t1
OC (collocation point) OC (collocation point)

OC (interpolation) OC (interpolation)
110 110
FD with 50 elements FD with 50 elements
rout = 0.00825 0.013813 rend = 0.021412 rout = 0.00825 0.013813 rend = 0.021412
r [m] r [m]
(c) (d)
Figure A.1: Temperature profiles for a step change in boundary value T (rout ) from 120 to
140 ◦ C generated by solving the one-dimensional nonlinear dynamic heat transfer problem in
radial direction of the PCM shell with orthogonal collocation (OC, solid line) or finite differ-
ences (FD, dashed line) discretizations.
In (a), with a single interior collocation point, the OC approximation is far off the FD bench-
mark trajectory even at the collocation point. By placing an additional collocation point at
rend , OC and FD approximations agree fairly well, and also the (geometric) mean tempera-
ture Tmean of the OC solution asymptotically converges to the FD benchmark.
Note that, for OC, the temperatures only at the collocation points fulfill the energy balance
equations, and continuous lines are obtained from inter-/extrapolating the trial functions in
eq. (13).
Figures (c) and (d) show the temperature profiles computed by OC with three and four col-
location points, respectively. Also see the error analysis in table A.1.
32
Appendix B. Transformation of differential term into standard form
Normalization of Ψ(r, x) in eq. (12) using eq. (11) and δ = rend − rout yields:

λP 1 ∂ ∂TP
Ψ(R, X) = 2 rout +δ(1−R)
δ rout +δ(1−R) ∂R ∂R

λP 1 ∂ 2 TP ∂ ∂TP
= 2 (rout +δ) −δ R
δ rout +δ(1−R) ∂R2 ∂R ∂R

λP 1 (rout +δ) ∂ ∂TP (rout +δ) ∂TP R ∂ ∂TP
= R − −δ R
(∗) δ 2 rout +δ(1−R) R ∂R ∂R R ∂R R ∂R ∂R

λP 1 1 ∂ ∂TP (rout +δ) ∂TP
= 2 (rout +δ−δR) R −
δ rout +δ(1−R) R ∂R ∂R R ∂R
" #
λP 1 ∂ ∂TP (rout +δ) 1 ∂TP
= 2 R −
δ R ∂R ∂R rout +δ(1−R) R ∂R
with the last line being the expression in eq. (15a).

∂ 2 TP 1 ∂ ∂TP 1 ∂TP
Note that in (∗) we have used = R − .
∂R2 R ∂R ∂R R ∂R
Appendix C. Derivation of the adapted matrix Ã(2)

In [38], the quadratic trial function evaluated at the collocation point j is
defined as:
N
X +2
yj := y(xj ) = x2i−2
j di (C.1)
i=1
which can be written for all j = 1, · · · , N + 2 collocation points as:

y = Qd (C.2)
using the notation y := (y1 , ..., yN +2 )T , d := (d1 , ..., dN +2 )T , and Q = (qji )ji with
qji := x2i−2
j .
For the first derivatives one obtains:
N
X +2 N
X +2
d d 2i−2
y(xj ) = xj di = (2i − 2) x2i−3
j di (C.3)
dx i=1
dx i=1
and, in the same way, for the expression in eq. (20c):

N
X +2 N
X +2
1 d 1 d 2i−2
zj := y(xj ) = xj di = (2i − 2) x2i−4
j di (C.4)
xj dx i=1
x j dx i=1
33
which can be written in matrix notation as
z = Cd = CQ−1 y = Ã(2) y (C.5)
with z = (z1 , ..., zN +2 )T and using d = Q−1 y from eq. (C.2). Then, the matrix
Ã(2) := CQ−1 is the adapted coefficient matrix used in eq. (20c).
Appendix D. Derivation of w(2) and the adapted vector w̃(2)
In Finlayson [41], the evaluation of integrals is based on the following quadra-

ture formula:
Z 1 N
X +2
a−1
f (x)x dx = wj f (xj ) (D.1)
0 j=1
(2) (2)
To determine wi and w̃i , evaluate eq. (D.1) for f (x) = x2i−2 using a = 2
and a = 1, respectively:
Z 1 N
X +2
1
x2i−2 xa−1 dx = wj x2i−2
j = =: fi , (D.2)
0 j=1
2i − 2 + a
wQ = f , w = fQ−1
using the notation f := (f1 , ..., fN +2 )T , w := (w1 , ..., wN +2 ), and Q = (qji )ji
with qji = x2i−2
j .
Appendix E. Estimation of non-measurable states
The accuracy of estimations of non-measurable PCM temperature states

is assessed by in silico studies. To this end ’artificial’ experimental data has
been generated by simulations with dMod for all states x, see eq. (27). To
test/benchmark the EKF in a more realistic setting with modeling errors differ-
ent model parameters in dMod as in the observer model rMod have been used.
However, the observer EKF-rMod is fed with measurable data only (meas).
The estimations of the measurable and the non-measurable states, generated by
EKF-rMod, are compared with the ’artificial’ experimental data. Note that the
measurable PCM temperature states in y correspond to the four temperature
points at the outer edge of the PCM shell: PCM 1, PCM 2, PCM 3, PCM 4.
The non-measurable states, noted as set-difference x\y, correspond to the four
temperature points inside the PCM shell: PCM 5, PCM 6, PCM 7, PCM 8.
See also table 1 for the notations.
The results are shown in figure E.2(a),(b). The first subfigures shows temper-
atures at the first axial position, at the outer edge (PCM 1) and inside (PCM 5)
of the PCM shell of the single tube. The second subfigures show temperatures
at the second axial position, at the outer edge (PCM 2) and inside (PCM 6) of
the PCM shell of the single tube. The third subfigures show the temperature
34
differences between PCM temperature states obtained from dMod and rMod
(figure E.2(a)) and dMod and EKF-rMod (figure E.2(b)).
While the simulations by the reduced model rMod show large deviations to
the measurements (figure E.2(a)), the EKF-rMod observer gives good estimates
for measurable and non-measurable states (figure E.2(b)). Analysis of all PCM
temperature states (PCM 1-8) reveals that the largest deviations always exist for
PCM 5 and PCM 1. For dMod and rMod in figure E.2(a) the maximal absolute
differences are 23.3 K and 17.3 K, respectively. Corresponding average absolute
errors are 0.03 K and 0.02 K. For dMod and EKF-rMod in figure E.2(b) the
maximal absolute differences are 5.1 K and 2.5 K, respectively. Corresponding
average absolute errors are 0.003 K and 0.007 K. For the studied scenario and
the assumed modeling errors it can be concluded that the observer produces
convincing results.
35
160
PCM # (dMod)
140 PCM 5 PCM # (rMod)
T [/ C]
120
100 PCM 1
80
160
PCM # (dMod)
140 PCM 6
PCM # (rMod)
T [/ C]
120
100 PCM 2
80
20
" PCM 1
10 " PCM 5
"T [K]
0
-10 " PCM 2
-20 " PCM 6
0 50 100 150 200 250 300

time [min]
(a)
160
PCM # (dMod)
140 PCM 5 PCM # (EKF-rMod)
T [/ C]
120
100 PCM 1
80
160
PCM # (dMod)
140 PCM 6
PCM # (EKF-rMod)
T [/ C]
120
100 PCM 2
80
20 " PCM 1
" PCM 5
10
"T [K]
0
-10 " PCM 2
" PCM 6
-20
0 50 100 150 200 250 300
time [min]
(b)
Figure E.2: (a) Predictions generated by simulations (dMod), (b) state estimations gen-
erated by the observer (EKF-rMod), and corresponding measurements (meas) for dynamic
experiment design (input) of HTF inlet temperature and a constant total mass flow of
ṁH,total = 1.09 kg/sec (same as in figure 9).
Measurements have been generated by simulations with dMod.
36
Bibliography
[1] A. d. Gracia and L. F. Cabeza. Phase change materials and thermal energy
storage for buildings. Energy and Buildings, 103:414–419, 2015.
[2] T. Nomura, N. Okinaka, and T. Akiyama. Technology of latent heat storage
for high temperature application: A review. ISIJ International, 50(9):1229–
1239, 2010. ISSN 0915-1559.
[3] A. Arteconi, N. J. Hewitt, and F. Polonara. State of the art of thermal
storage for demand-side management. Applied Energy, 93:371–389, 2012.
[4] M. A. Sharif, A. A. Al-Abidi, S. Mat, K. Sopian, M. H. Ruslan, M. Y.
Sulaiman, and M. A. M. Rosli. Review of the application of phase change
material for heating and domestic hot water systems. Renewable and Sus-
tainable Energy Reviews, 42:557–568, 2015.
[5] M. Mastani Joybari, F. Haghighat, J. Moffat, and P. Sra. Heat and cold
storage using phase change materials in domestic refrigeration systems:
The state-of-the-art review. Energy and Buildings, 106:111–124, 2015.
[6] B. Cárdenas and N. León. High temperature latent heat thermal energy
storage: Phase change materials, design considerations and performance
enhancement techniques. Renewable and Sustainable Energy Reviews, 27:
724–737, 2013. ISSN 13640321.
[7] Z. Youssef, A. Delahaye, L. Huang, F. Trinquet, L. Fournaison, C. Poller-
berg, and C. Doetsch. State of the art on phase change material slurries.
Energy Conversion and Management, 65:120–132, 2013.
[8] M. Liu, W. Saman, and F. Bruno. Review on storage materials and thermal
performance enhancement techniques for high temperature phase change
thermal storage systems. Renewable and Sustainable Energy Reviews, 16
(4):2118–2132, 2012.
[9] W. J. Cole, K. M. Powell, and T. F. Edgar. Optimization and control of
thermal energy storage systems. Reviews in Chemical Engineering, 28(2-3):
81–99, 2012.
[10] K. M. Powell, W. J. Cole, U. F. Ekarika, and T. F. Edgar. Optimal chiller
loading in a district cooling system with thermal energy storage. Energy,
50:445–453, 2013.
[11] T. F. Edgar and K. M. Powell. Energy intensification using thermal storage.
Current Opinion in Chemical Engineering, 9:83–88, 2015.
[12] A. Castell and C. Solé. Design of latent heat storage systems using phase
change materials (PCMs). In L. F. Cabeza, editor, Advances in thermal en-
ergy storage systems: methods and applications, pages 307–324. Woodhead
Publishing, 2014.
37
[13] E. Oró, A. d. Gracia, A. Castell, M. M. Farid, and L. F. Cabeza. Re-
view on phase change materials (PCMs) for cold thermal energy storage
applications. Applied Energy, 99:513–533, 2012.
[14] S. Tiari, S. Qiu, and M. Mahdavi. Discharging process of a finned heat
pipe–assisted thermal energy storage system with high temperature phase
change material. Energy Conversion and Management, 118:426–437, 2016.
[15] H. J. Strumpf, V. Avanessian, and R. Ghafourian. Design analysis and
containment canister life prediction for a brayton engine solar receiver for
space station. Journal of Solar Energy Engineering, 116(3):142–147, 1994.
[16] F. Hengstberger, C. Zauner, K. Resch, S. Holper, and M. Grobbauer. High
temperature phase change materials for the overheating protection of facade
integrated solar thermal collectors. Energy and Buildings, 124:1–6, 2016.
[17] T. Barz, C. Zauner, D. Lager, D. C. López C., F. Hengstberger, M. N.
Cruz Bournazou, and K. Marx. Experimental analysis and numerical mod-
eling of a shell and tube heat storage unit with phase change materials.
Industrial & Engineering Chemistry Research, 55(29):8154–8164, 2016.
[18] S. D. Vaca Jiménez. Model to determine the state of charge of water as
a phase change material integrated in a household refrigerator. Master’s
thesis, Loughborough University, under the supervision of P. Eames and
W. S. Ranshuysen, Leicestershire, United Kingdom, 02.12.2013.
[19] A. Cuneo, M. L. Ferrari, M. Pascenti, and A. Traverso. State of charge
estimation of thermal storages for distributed generation systems. Energy
Procedia, 61:254–257, 2014.
[20] C. Hall, E. Glakpe, J. Cannon, and T. Kerslake. Thermal state-of-charge
in solar heat receivers. In 36th Aerospace Sciences Meeting & Exhibit, 1998.
doi: 10.2514/6.1998-1017.
[21] G. P. Henze, D. E. Kalz, S. Liu, and C. Felsmann. Experimental analysis
of model-based predictive optimal control for active and passive building
thermal storage inventory. HVAC&R Research, 11(2):189–213, 2005.
[22] M. S. Naghavi, K. S. Ong, M. Mehrali, I. A. Badruddin, and H. S. C.
Metselaar. A state-of-the-art review on hybrid heat pipe latent heat storage
systems. Energy Conversion and Management, 105:1178–1204, 2015.
[23] G. Steinmaurer, M. Krupa, and P. Kefer. Development of sensors for mea-
suring the enthalpy of PCM storage systems. Energy Procedia, 48:440–446,
2014.
[24] T. Bauer, D. Laing, and W.-D. Steinmann. Feasibility of the new PCM
measurement system based on the ’electric resistance approach’ in lab scale
demonstrated. Deliverable report 15.4 - Solar Facilities for the European
Research Area (SFERA) Project, German Aerospace Center (DLR), Ger-
many, 2012.
38
[25] T. Kreuzinger, M. Bitzer, and W. Marquardt. State estimation of a strat-
ified storage tank. Control Engineering Practice, 16(3):308–320, 2008.
[26] Q. Mao. Recent developments in geometrical configurations of thermal en-
ergy storage for concentrating solar power plant. Renewable and Sustainable
Energy Reviews, 59:320–327, 2016.
[27] F. Agyenim, N. Hewitt, P. Eames, and M. Smyth. A review of materials,
heat transfer and phase change problem formulation for latent heat ther-
mal energy storage systems (LHTESS). Renewable and Sustainable Energy
Reviews, 14(2):615–628, 2010.
[28] C. Zauner, F. Hengstberger, M. Etzel, D. Lange, R. Hofmann, and H. Wal-
ter. Experimental characterization and simulation of a fin-tube latent heat
storage using high density polyethylene as PCM. Applied Energy, 179:
237–246, 2016.
[29] A. Trp. An experimental and numerical investigation of heat transfer during
technical grade paraffin melting and solidification in a shell-and-tube latent
thermal energy storage unit. Solar Energy, 79(6):648–660, 2005.
[30] V. R. Voller, C. R. Swaminathan, and B. G. Thomas. Fixed grid techniques
for phase change problems: a review. International Journal for Numerical
Methods in Engineering, 30(4):875–898, 1990.
[31] H. Hu and S. A. Argyropoulos. Mathematical modelling of solidification
and melting: a review. Modelling and Simulation in Materials Science and
Engineering, 4(4):371–396, 1996.
[32] M. Muhieddine, E. Canot, and R. March. Various approaches for solving
problems in heat conduction with phase change. International Journal on
Finite Volumes, 6(1):1–20, 2009.
[33] G. Ziskind. Modelling of heat transfer in phase change materials (PCMs)
for thermal energy storage systems. In L. F. Cabeza, editor, Advances in
thermal energy storage systems: methods and applications, pages 307–324.
Woodhead Publishing, 2014.
[34] A. Samarskii, P. N. Vabishchevich, O. P. Iliev, and A. G. Churbanov. Nu-
merical simulation of convection/diffusion phase change problems - a re-
view. International Journal of Heat and Mass Transfer, 36(17):4095–4106,
1993.
[35] U. Gaur and B. Wunderlich. Heat capacity and other thermodynamic
properties of linear macromolecules. ii. polyethylene. Journal of Physical
and Chemical Reference Data, 10(1):119, 1981.
[36] A. Y. Uzan, Y. Kozak, Y. Korin, I. Harary, H. Mehling, and G. Ziskind. A
novel multi-dimensional model for solidification process with supercooling.
International Journal of Heat and Mass Transfer, 106:91–102, 2017.
39
[37] NIST/SEMATECH. e-handbook of statistical methods, 1 2015. URL http:
//www.itl.nist.gov/div898/handbook/.
[38] B. A. Finlayson. Nonlinear Analysis in Chemical Engineering. McGraw-
Hill, New York, 1980.
[39] R. G. Rice and D. D. Do. Applied mathematics and modeling for chemical
engineers. John Wiley & Sons, 2012. ISBN 1118343026.
[40] S. Quoilin, A. Desideri, J. Wronski, I. Bell, and V. Lemort. Thermocycle:
A modelica library for the simulation of thermodynamic systems. In The
10th International Modelica Conference, March 10-12, 2014, Lund, Sweden,
Linköping Electronic Conference Proceedings, pages 683–692. Linköping
University Electronic Press, 2014.
[41] B. A. Finlayson. The method of weighted residuals and variational princi-
ples. ACADEMIC PRESS, INC., New York, 1972. ISBN 1611973236.
[42] N. L. Segall, J. F. MacGregor, and J. D. Wright. Collocation methods for
solving packed bed reactor models with radial gradients. The Canadian
Journal of Chemical Engineering, 62(6):808–817, 1984.
[43] J. B. Rawlings and B. R. Bakshi. Particle filtering and moving horizon es-
timation. Computers & Chemical Engineering, 30(10-12):1529–1541, 2006.
[44] A. Gelb. Applied Optimal Estimation. MIT Press, Cambridge, MA, 1974.
[45] R. Schneider and C. Georgakis. How to NOT make the extended Kalman
filter fail. Industrial & Engineering Chemistry Research, 52(9):3354–3362,
2013.
[46] R. E. Murray and D. Groulx. Experimental study of the phase change
and energy characteristics inside a cylindrical latent heat energy storage
system: Part 1 consecutive charging and discharging. Renewable Energy,
62:571–581, 2014.
[47] P. Kühl, M. Diehl, T. Kraus, J. P. Schlöder, and H. G. Bock. A real-time
algorithm for moving horizon state and parameter estimation. Computers
& Chemical Engineering, 35(1):71 – 83, 2011. ISSN 0098-1354.
[48] M. Dannemand, J. M. Schultz, J. B. Johansen, and S. Furbo. Long term
thermal energy storage with stable supercooled sodium acetate trihydrate.
Applied Thermal Engineering, 91:671–678, 2015.
40
View publication stats

2 D Model

Uploaded by

Copyright:

Available Formats

2 D Model

Uploaded by

Document Information

Original Description:

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

2 D Model

Uploaded by

Copyright:

Available Formats

See discussions, stats, and author profiles for this publication at: https://www.researchgate.

State and state of charge estimation for a latent heat storage

Article in Control Engineering Practice · March 2018

Tilman Barz Dominik Seliger

SEE PROFILE SEE PROFILE

Klemens Marx Andreas Sommer

SEE PROFILE SEE PROFILE

WHole Battery - Wärmeoptimierte Hochleistungs-Batterie View project

AdaMo Project View project

The user has requested enhancement of the downloaded file.

Tilman Barz1 , Dominik Seliger, Klemens Marx

Andreas Sommer, Sebastian F. Walter, Hans Georg Bock

1 Corresponding author: tilman.barz@ait.ac.at

Thermal energy storages (TES) have a great potential in energy conserva-

State of charge (SOC) of latent thermal storages: Undoubtedly, monitoring of

Measurement based approaches for estimating SOC: Reports on measurement

2. Shell and tube 2D LHTES Model

2.1. Container/storage configuration

The mathematical modeling of the conductive heat transfer in the PCM

Balance equation for heat transfer fluid (subscript H)

0 ∂TW (t, x) ∂TW (t, x)

Balance equation for phase change material (subscript P)

∂TP (t, r, x) ∂TP (t, r, x) ∂TP (t, r, x)

q̇H (t, x) = α(t, x) (TH (t, x) − TW (t, x)) (4a)

with α being the heat transfer coefficient function.

λP (t, r, x) r=rout ∂TP (t, r, x)

where ξ = ml /(ms + ml ) is the liquid weight (mass) fraction, or phase fraction,

3. Spatial discretization of the reduced single tube model (rMod)

3.1. Normalization of the independent space variables

Before normalization of space variables, the diffusion terms in x in eq. (2a)

3.2. Collocation scheme for the PCM

NR = 1 and NX = 4 define the number of interior collocation points as well

4. Details on the numerical implementation (rMod)

4.1. Normalized equations

The boundary condition for rMod in eq. (6) reads

λP (t, R, X) R=1 ∂TP (t, R, X)

4.2. Collocation applied to the PCM shell

4.3. Evaluation of integrals

5. Comparison of detailed and reduced single tube models dMod and

dMod represents a well accepted model complexity for numerical studies of

notation variable symbol

rMod yields a system of 12 continuous-time nonlinear state equations and 5

ẋ(t) = f (x(t), u(t), p) (26)

x = [TH (X1 ), TH (X2 ), TH (X3 ), TH (X4 ), (27)

160 HTF in (input)

160 PCM 1-4 (dMod)

120 PCM 1, 2,3,4

1 HTF flow (input)

145 rend =0.021 rend =0.021

140 0.014 135 0.014

6. The state observer

An overview of currently available observer methods for state estimation of

6.1. State observer design and observability

6.2. Numerical implementation

6.3. Experimental results

• rMod: predictions generated by simulations with rMod,

HTF out (rMod)

160 HTF in (input)