Renewable Energy 150 (2020) 891e906

Contents lists available at ScienceDirect

Renewable Energy
journal homepage: www.elsevier.com/locate/renene

Energy modeling of solar water heating systems with on-off control

and thermally stratified storage using a fast computation algorithm

nio Araújo*, Rui Silva
Anto
~o, Portugal
Faculdade de Engenharias e Tecnologias, Universidade Lusíada Norte, 4760-108, Vila Nova de Famalica

a r t i c l e i n f o a b s t r a c t

Article history: This work presents a simplified model for the rapid computation of the yearly solar fraction of direct
Received 8 September 2019 solar water heating systems using on-off control. Thermal stratification was included using a simple one-
Received in revised form dimensional multi-node model. A time-step dependency analysis showed that a time step of 0:05 h is a
6 December 2019
good compromise between accuracy and computation speed. The solar fraction increases with collector
Accepted 6 January 2020
flow rate when the flow rate is low. In fully-mixed storage, the solar fraction keeps increasing with flow
Available online 13 January 2020
rate, although with a decreasing rate of increase. However, in stratified storage, the solar fraction reaches
a maximum at an optimum flow rate, before it starts decreasing with flow rate. When the number of tank
keywords:
Solar water heating
nodes increases from 1 to 4, the maximum solar fraction increases 5e28 %; this increase is superior for
On-off control less efficient collectors and lower collector areas. In low-stratified systems, the optimum flow rate is the
Thermal stratification maximum allowed by the system. However, in stratified systems, the optimum flow rate is reduced to
Simulation values of 0.006e0:016 m3 h
1 per square meter of collector area. Unless the tank walls are covered by a
rather thick layer of thermal insulation (about 0:2 m), storage tank losses cannot be ignored.
© 2020 Elsevier Ltd. All rights reserved.

1. Introduction continuously varied so as to maintain at a preset fixed value the

temperature of the fluid leaving the collector.
The importance of domestic solar water heating systems may be Active solar water heating may be further classified into direct
appreciated by the fact that, in 2014, small-scale systems in single- systems, if the solar collector and the storage tank share the same
family houses contributed to 68 % of the solar thermal energy working fluid, or indirect systems, if the collector fluid is different
produced worldwide [1]. from that of the storage tank. In direct systems, heat is transferred
Solar water heating systems comprise the following main ele- to the storage tank by direct fluid mixing; in indirect systems, heat
ments: a collection unit that converts solar radiation into useable is transferred to the storage tank by means of a heat exchanger [3].
thermal energy, a storage tank for the accumulation of thermal Solar water heating systems are dynamic and complex; their
energy, and an auxiliary heater to compensate for energy shortfalls. performance is dependent on time-immutable parameters, such as
These elements are interconnected by insulated pipelines, through the physical and geometrical characteristics of system equipment,
which a working fluid flows. and on time-variable parameters, such as the weather conditions
Solar water heating systems are classified into active systems, and consumption load. As a result, the modeling of solar water
when a pump is used to force the working fluid through the solar heating systems involves the solution of various interrelated time-
collection unit, or passive systems, when fluid flow occurs by nat- dependent differential equations, which are characterized by non-
ural convection [2,3]. In active solar heating systems, two basic linearity and transient behavior. Several time-marching numerical
schemes are used to control fluid flow through the solar collector techniques, often involving numerous iterations per time step, have
[2,4]: on-off control, where the flow rate is fixed and either turned been used to solve the differential equations that model different
on or off, depending on the temperature of the fluid exiting the solar water heating systems [5e10]. Third-party computer pack-
solar collector; proportional control, where the flow rate is ages, such as Simulink from MathWorks [11] and the general-
purpose simulation program TRNSYS [12], have also been
frequently used to aid the modeling of solar water heating systems
[13e26].
* Corresponding author.
As simulations based on a specific weather database may
E-mail address: antonio.araujo@hotmail.com (A. Araújo).

https://doi.org/10.1016/j.renene.2020.01.026
0960-1481/© 2020 Elsevier Ltd. All rights reserved.
892 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906

Nomenclature AS;k Storage area of node k (m2 )

C_ B Consumption/storage heat capacity flow rate
b Collector tilt (rad) (kW + C
1 )
g Collector azimuth (rad) CC Daily consumption heat capacity (kWh + C
1 )
d Solar declination (rad) CS;k Storage heat capacity of node k (kWh + C
1 )
dF Solar fraction error DS Storage diameter (m)
dF max Relative decrease of maximum solar fraction F Yearly solar fraction
dT ref
S Maximum relative difference between storage FU Heat loss coefficient of collector (kW m
2 + C
1 )
temperatures on first and second simulation years F max Maximum yearly solar fraction
dTS;k Relative difference between storage temperatures of F max;0 Maximum solar fraction of perfectly insulated
node k on first and second simulation years storage
h Collector efficiency F ref Reference solar fraction
h0 Optical efficiency of collector G Solar radiation energy (kWh m
2 )
q Empirical parameter of mains temperature model G_ Solar radiation power (kW m
2 )
(rad) G0;z Horizontal extraterrestrial solar radiation energy
r Density of working fluid (kg m
3 ) (kWh m
2 )
rg Reflectivity of the surrounding surfaces Gb;z Horizontal beam solar radiation energy (kWh m
2 )
f Latitude (rad) Gd;z Horizontal diffuse solar radiation energy (kWh m
2 )
u Solar angle (rad) N Number of simulation days
u0 Initial solar angle (rad) QA Consumption of auxiliary energy (kWh)
u1 Final solar angle (rad) Rb Ratio of collector to horizontal beam radiation
usr Solar angle at sunrise (rad) Sb Collector beam radiation factor
uss Solar angle at sunset (rad) Sb;z Horizontal beam radiation factor
dt Time step (h) TB Consumption/storage temperature (+ C)
dt ref Reference time step (h) TC Consumption temperature (+ C)
dxS Storage wall thickness (m) TE Outdoor temperarure (+ C)
dTE Yearly outdoor temperature amplitude (+ C)
TE Yearly mean outdoor temperature (+ C)
dTP Collector temperature difference (+ C)
TI Collector inlet temperature (+ C)
SQA Yearly consumption of auxiliary energy (kWh)
TM Mains temperature (+ C)
SQC Yearly energy consumption (kWh)
TN Indoor temperarure (+ C)
a Anisotropy index
TP Production temperature (+ C)
b Number of storage nodes
T max
P Maximum production temperature (+ C)
c Specific heat capacity of working fluid
TS;k Storage temperature of node k (+ C)
(kWh kg
1  C
1 )
T 0S;k Initial storage temperature of node k (+ C)
f Modulating factor
T 1S;k Final storage temperature of node k (+ C)
hS Storage external coefficient of convection
T new Storage temperature of node k on second simulation
(kW m
2 + C
1 ) S;k
year (+ C)
k Storage node
T old Storage temperature of node k on first simulation
kS Storage wall thermal conductivity (kW m
1 + C
1 ) S;k
year (+ C)
lS Storage length (m)
TU Temperature of storage location (+ C)
n Day of the year
US Storage heat transfer coefficient (kWh m
2 + C
1 )
p Production node
pC Number of consumers V_ B Consumption/storage volume flow rate (m3 h
1 )
r1 Empirical parameter of mains temperature model V_ P Production volume flow rate (m3 h
1 )

r2 Empirical parameter of mains temperature model V_ P Optimum production volume flow rate (m3 h
1 )
rS Storage length-to-diameter ratio max
V_ P Maximum production volume flow rate (m3 h
1 )
t Solar time (h)
VS Storage volume (m3 )
t0 Initial solar time (h)
VS;k Storage volume of node k (m3 )
t1 Final solar time (h)
SG Single-glazed
vC Daily consumption volume (m3 )
UG Unglazed
xC Hourly consumption fraction
A Collector area (m2 )

provide quite different results from those based on real weather above a certain level of system detail becomes redundant due the
data, and as simulation results based on a predefined consumption difference between real and simulated values of weather data and
profile may be quite different from those that would be obtained consumption profile. On the other hand, due to the inherent
from real consumption data, Araújo and Pereira [27,28] concluded complexity of the solar water heating systems, added system de-
that, in general, it is not worth to build overly complex models to tails requires highly sophisticated numerical methods, resulting in
simulate the long-term performance of solar water heating systems very time-consuming iterative computations. As a result, Araújo
due to the unavoidable uncertainty introduced by the weather data and Pereira developed simplified solar water heating models,
and consumption profile. This is so because the added accuracy that which promote linearity and exclude iterative processes and im-
would be introduced by the fine-tuning of system parameters plicit solutions, for systems using both on-off [27] and proportional
 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906 893

[28] control schemes. The simplicity of these model enables the parameters (e.g., fluid inlet flow rate and temperature); modeling
rapid computation of the long-term performance of solar heating methods for thermal stratification that predict the transient
systems, making these models ideal for optimization procedures, in behavior of the temperature distribution within the storage tank in
which an objective cost function has to be computed several times one, two, or three dimensions.
before an optimal solution is reached. However, while the added In direct systems with thermal stratification, the most impor-
complexity of fine system-detail modeling may turn out to be tant operational parameter is likely to be the inlet/outlet fluid flow
redundant if only the absolute long-term performance of real sys- rate between the collection unit and the storage tank, as high flow
tems is to be predicted, the effort put into the modeling of system rates tend to cause the destruction of thermal stratification due to
details may become very important in comparative studies aiming excessive fluid mixing [31,39]. Collector flow rates of around
to evaluate the relative effect of varying system parameters on 0:05 m3 h
1 per square meter of useful collector area are not un-
system performance. In fact, it is advantageous if a constant common if thermal stratification is not taken into account. How-
weather data set and a constant consumption profile are shared ever, optimal values of about one fifth of this value (i.e., about
among different simulations of comparative studies. 0:01 m3 h
1 ) are typically required in order to maintain high de-
One aspect that was not taken into account by Araújo and Per- grees of stratification [2,41,42].
eira [27,28] was the phenomenon of thermal stratification of the Stratified tanks are very difficult to model, except in the limit of
fluid within the storage tank, i.e., the storage fluid was assumed to perfect stratification. However, as thermal stratification in real
be perfectly mixed. Basically, the concept of thermal stratification tanks is always between a perfectly mixed and a perfectly stratified
lies in the fact that, as the density of liquids (e.g., water) increases tank, in many practical problems, it may not be necessary more
with increasing temperature, warmer parts of the working fluid in than these two stratification limits to characterize the degree of
the storage tank tend to move upwards relative to colder parts of stratification in real tanks [29,39]. Moreover, even though complex
the fluid, i.e., an increasing temperature gradient is established models in two or three dimensions may take into account more
from the bottom to the top of the fluid within the storage tank factors affecting the performance of thermally stratified systems,
[2,29,30]. Thus, it is assumed that any incoming fluid will be simpler one-dimensional models may be advantageous because
directed to the level at which its density matches the density of the they are computationally more efficient [45]. Consequently, a
surrounding fluid [29e31]. Consequently, in order to promote highly complex model of temperature stratification within a tank is
stratification, hot fluid inlets are normally located at the top of the often not necessary for predicting system performance [39].
tank, whereas cold fluid inlets are located at the bottom of the The objective of the present work is to extend the fast compu-
storage tank. tation model developed by Araújo and Pereira [27] to include
Thermal stratification increases the efficiency of solar water thermal stratification in the storage tank for systems with direct
heating systems due to two main reasons [29,31e33]: one the one fluid mixing. The purpose of the model is to provide a fast simu-
hand, by placing the outlet to the solar collector at the bottom of the lation means of computing a reasonably accurate estimate of the
storage tank, the efficiency of the solar collector increases due to long-term performance of direct solar water heating systems using
the lower temperature at the collector inlet; on the other hand, by the on-off control scheme. As discussed earlier, if only the absolute
placing the outlet to the consumption circuit at the top of the performance of solar water heating systems is required, there is no
storage tank, the demand for auxiliary energy is reduced due to the reason to build overly complex models because they tend to be
higher temperature of the fluid available for consumption. computationally inefficient. Consequently, the thermal stratifica-
By a comparison between fully mixed and stratified tanks used tion model should be easily integrated with the other parts of the
in solar water heating, system efficiency rises in the range of solar thermal model and simple and fast to compute.
5e20 % or higher have been reported in the literature: Sharp and In the following section (Section 2), the solar water heating
Loehrke [34] reported improvements in system performance due to model is formulated: physical models are developed for a stratified
thermal stratification in the range of 5e15 %; van Koppen et al. [35] storage tank, solar collection unit, and consumption circuit; then,
stated that an increase in heat gain of 5e10 % is likely to be ex- these three physical models are integrated into a numerical pro-
pected; Veltkanp [36] showed that the proper utilization of ther- cedure to compute long-term system performance; solar radiation
mally stratified storage may enhance the output of solar thermal and mains water temperature models are also formulated. Model
systems between 10 and 20 %; performance improvements of simulation parameters are defined in Section 3, which include
5e20 % were reported by Cole and Belinger [37]; Wuestling et al. weather data, design parameters, a consumption profile, system
[38] reported performance improvements in the range of 12e15 %; performance, and numerical parameters. The numerical procedure
Duffie [39] stated that system performance improvements in the is validated by means of a time-step dependency analysis and the
range of 10e20 % may be obtained in highly stratified systems; study of a stopping criterion to guarantee periodic convergence. In
Ghaddar [40] reported an increase of up to 20 % in the energy Section 4, a sensitivity analysis is performed for the following pa-
delivered when stratification is employed in the storage tank; rameters: flow rate through the collector, degree of stratification in
Hollands and Lightstone [41] reported the highest improvements the storage tank, solar collector area, storage tank volume, and tank
due to thermal stratification, reporting efficiency increases as high thermal losses. Through the sensitivity analysis, the solar heating
as 38 %. However, almost all studies quantifying the performance model can be evaluated against published results. Finally, in Section
gains achieved by solar water heating systems with thermally 5, model results are summarized and briefly discussed.
stratified storage over fully mixed systems are over 30 years old
[42]. 2. Solar water heating model
According to the main state-of-the-art papers from the last 10
years [30,31,43,44], most published material on thermal stratifi- The main objective of the solar thermal model developed herein
cation has focused mainly on three topics: methods to quantify the is the computation of the yearly consumption of auxiliary energy
degree of stratification, which is typically characterized by a single (SQA , kWh) that has to be supplied due to insufficient collection
dimension or dimensionless parameter; stratification improve- and/or storage of solar energy in active solar water heating systems.
ment parameters, typically divided into physical/geometrical pa- The computation process runs through successive discrete time
rameters (e.g., tank size and geometry, placement and geometry of steps (dt, h), and, at each time step, the consumption of auxiliary
inlet ports, and thermal insulation characteristics) and operational energy (QA , kWh) is computed, so that the values of QA from all
894 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906

Fig. 1. Solar thermal system.

time steps can be added up, resulting in yearly consumption SQA .

As indicated above and shown in Fig. 1, a solar water heating
system is formed by three mains parts: the production circuit, the
storage tank, and the consumption circuit. In the present work, it is
assumed that the working fluid of the production circuit is the same
as that within the storage tank and that flowing in the consumption
circuit. In the production circuit, the solar collector absorbs solar
radiation and transfers it to the working fluid. The on-off scheme is
used to control fluid flow through the production circuit. The
thermal energy of the production fluid is subsequently transferred
to the storage tank by direct fluid mixing, and the energy accu-
mulated in the storage fluid is transferred to the consumption cir-
cuit whenever it is required. When the temperature of the
consumption fluid is not sufficiently high, the fluid has to be heated
up by means of an auxiliary heater.
Sections 2.1e2.3 present the development of the physical
models for the storage, production, and consumption systems,
respectively. Sections 2.5 and 2.6 present a physical model for the
useful solar radiation and an empirical model for the mains water
temperature, respectively. Section 2.4 presents the numerical
implementation of the physical models. The numerical model was
implemented using Julia, a high-level general-purpose dynamic
programming language for high-performance numerical analysis
and computational science [46,47].
Fig. 2. Energetic fluxes in stratified storage tank.

2.1. Storage system

the following energy balances can be formulated: if 1  k  p
1,
The modeling of the storage system was based on the multi-
node approach presented by Duffie and Beckman [2] for ther-
mally stratified tanks. As shown in Fig. 2, the tank is divided into b   dTS;k
C_ B TS;kþ1
C_ B TS;k
US AS;k TS;k
TU ¼ CS;k ; (1a)
vertically stacked sections (nodes) of volume VS;k (m3 ), where k ¼ dt
1; 2; …; b, so that the total volume of the storage tank
if k ¼ p and psb,
X
b
VS ¼ VS;k :  
k¼1
C_ P TP þ C_ B TS;kþ1
C_ P TS;k
C_ B TS;k
US AS;k TS;k
TU
dTS;k
The fluid within each volume VS;k is fully mixed at temperature ¼ CS;k ; (1b)
TS;k (+ C). dt
At the bottom (k ¼ b) of the storage tank, the storage fluid flows
if p þ 1  k  b
1,
to the production circuit at temperature TS;b and heat capacity flow
rate C_ P ¼ rV_ P c, where r (kg m
3 ) and c (kWh kg
1  C
1 ) are the  
density and specific heat capacity of the working fluid, respectively, C_ P TS;k
1 þ C_ B TS;kþ1
C_ P TS;k
C_ B TS;k
US AS;k TS;k
TU
and V_ P (m3 h
1 ) is the volume flow rate in the production circuit. dTS;k
The consumption fluid enters the bottom of the storage tank at ¼ CS;k ; (1c)
dt
mains temperature TM (+ C) and heat capacity flow rate C_ B ¼ rV_ B c,
where V_ B (m3 h
1 ) is the volume flow rate of the consumption fluid if k ¼ b and psb,
entering the storage tank. At the top (k ¼ 1) of the storage tank, the
storage fluid flows to the consumption circuit at temperature TS;1  
and heat capacity flow rate C_ B . The production fluid may enter the C_ P TS;k
1 þ C_ B TM
C_ P TS;k
C_ B TS;k
US AS;k TS;k
TU
storage tank at any node p, where 1  p  b, at production tem- dTS;k
perature TP (+ C) and heat capacity flow rate C_ P . ¼ CS;k ; (1d)
dt
The energy rates in and out of node k depend on the relative
position of k with respect to p and of p with respect to b (Fig. 2), and if k ¼ p ¼ b,
 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906 895

heated from temperature TB (+ C) up to the temperature required for

  dTS;k consumption (TC , + C) when TB < TC , but it must be mixed with
C_ P TP þ C_ B TM
C_ P TS;k
C_ B TS;k
US AS;k TS;k
TU ¼ CS;k ;
dt colder consumption fluid at temperature TM when TB > TC .
(1e) When TB < T , flow rate C_ B ¼ C_ , where C_ ¼ rV_ c, and V_
C C C C C
(m3 h
1 ) is the volume flow rate required for consumption. The
where CS;k ¼ rVS;k c is the storage heat capacity, and t (h) represents energy supplied by the auxiliary heater
time. Term US AS;k ðTS;k
TU Þ represents the heat loss though the
tank walls at node k, where US (kWh m
2 + C
1 ) is the global heat
QA ¼ C_ C ðTC
TB Þdt: (6)
transfer coefficient of the tank walls, AS;k (m2 ) is the heat transfer
area of node k, and TU (+ C) is the temperature outside the storage When TB  TC , auxiliary energy QA ¼ 0 kWh. A mass and an
tank. The tank is assumed to be located in an interior uncondi- energy balance to the fluid mixing process results in the following
tioned space, whose temperature may be approximated to the expression for the flow rate of the fluid entering the storage tank:
mean of the indoor conditioned temperature (TN , + C) and the
outdoor temperature (TE , + C) [48]: T
TM _
C_ B ¼ C C : (7)
TB
TM C
T þ TE
TU ¼ N : (2) As mention in Section 2.1, the storage fluid flows directly from
2
the top of the storage tank to the consumption circuit, so that TB ¼
TS;1 .

2.2. Production circuit

As shown in Fig. 1, the production fluid enters the solar collector

2.4. Numerical procedure
at volume flow rate V_ P and temperature TI (+ C), leaving it at a
higher temperature TP due to the solar radiation power flux (G, _
For the differential equations (1a)e(1e), approximate numerical
kW m
2 ) reaching the collector, so that the following energy bal-
solutions were developed for the set of discrete time instances
ance can be established:
separated by time steps dt. The forward Euler method [51] was
_ ¼ C_ ðT
T Þ;
hGA (3) applied to calculate approximate storage temperatures at the end of
P P I
each time step, so that, for time step dt, TS;k zT 0S;k , and
where h is the efficiency of the collector, and A (m2 ) is the useful
1 0
collector area. dTS;k T S;k
T S;k
For simplicity, a linear relation is used to model collector effi- z ;
dt dt
ciency [2,49,50]:
where T 0S;k and T 1S;k are the values of the storage temperatures at the
dTP beginning and end of dt, respectively, of nodes k ¼ 1; 2; …; b,
h ¼ h0
FU _ ; (4)
G resulting in the following explicit solutions for T 1S;k : if 1  k  p
1,

where vertical intercept h0 is the optical efficiency, and slope FU    

dt
(kW m
2 + C
1 ) is the heat loss coefficient of the collector. In the
1
TS;k ¼ C_ B TS;kþ1
0
þ US AS;k TU
C_ B þ US AS;k TS;k
0 0
þ TS;k ;
CS;k
present work, temperature difference dTP follows the European
practice [2,49], i.e., it is the difference between the mid-range (8a)
temperature of the fluid in the collector and the outdoor
if k ¼ p and psb,
temperature:
 
TP þ TI dt
dTP ¼
TE :
1
TS;k ¼ C_ P TP þ C_ B TS;kþ1
0
þ US AS;k TU
C_ P þ C_ B
2 CS;k
 
In on-off control, flow rate C_ P is fixed, and fluid flow is either 0 0
þ US AS;k TS;k þ TS;k ; (8b)
turned on or off, depending on the value of temperature TP , which
must be, therefore, continuously monitored. Temperature TP can be
computed by simultaneously solving Eqs. (3) and (4) with respect if p þ 1  k  b
1,
to TP :  
dt
1
TS;k ¼ C_ P TS;k
1
0
þ C_ B TS;kþ1
0
þ US AS;k TU
C_ P þ C_ B
h G_ þ FUðTE
TI Þ CS;k
TP ¼ 2A 0 þ TI : (5)  
2C_ þ FUA
P 0 0
þ US AS;k TS;k þ TS;k ; (8c)
The flow of the production fluid is turned on whenever TI < TP <
T max
P and turned off otherwise, where T max
P (+ C) is the maximum
temperature of the production fluid. if k ¼ b and psb,
As mention in Section 2.1, the storage fluid flows directly from  
the bottom of the storage tank to the production circuit, so that TI ¼ dt
1
TS;k ¼ C_ P TS;k
1
0
þ C_ B TM þ US AS;k TU
C_ P þ C_ B
TS;b . CS;k
 
0 0
þ US AS;k TS;k þ TS;k ; (8d)
2.3. Consumption circuit

As shown in Fig. 1, the fluid of the consumption circuit must be if k ¼ p ¼ b,

896 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906

   
dt
1
TS;k ¼ C_ P TP þ C_ B TM þ US AS;k TU
C_ P þ C_ B þ US AS;k TS;k
0
CS;k
0
þ TS;k :
(8e)
The computation procedure used to calculate storage tempera-
tures T 1S;1 ; T 1S;2 ; …; T 1S;b is shown in Fig. 3.
The forward Euler method was employed even though its
convergence is restricted to short time steps, since it produces
rather simple close-form solutions, as stated by Eqs. (8a)e(8e).
Although the backward Euler method is absolutely stable (i.e., it
converges for any time step), this approach was disregarded
because it yields overly complex implicit solutions for T 1S;k, whose
complexity increases with increasing b, that must be found by some
root-finding algorithm [51]. A study of the effect of the size of dt on
the solutions obtained using the forward Euler method is presented
in Section 3.5.1.
For the computation of Eqs. (8a)e(8e), three quantities must be
computed in advance: production temperature TP , fluid entrance
node p, and flow rate C_ B . The values of TP and p depend on the
conditions in the production circuit and in the storage tank; the
value of C_ B depends on the conditions in the consumption circuit
and in the storage tank.
Following the forward Euler method [51], temperatures TI and
TB must be equated to storage temperatures TS;b and TS;1 , respec-
tively, at the beginning of time step dt, i.e., TI ¼ T 0S;b , and TB ¼ T 0S;1 .
In thermally stratified systems, the production fluid tends to
Fig. 3. Flow chart of the computation procedure of storage temperatures T 1S;1 ; T 1S;2 ; …; enter the storage tank at the level at which the temperatures of the
T 1S;b .
production and storage fluids are similar [2,29]. Therefore, at each
time step dt, entrance node p is the topmost node for which pro-
duction temperature TP is higher than the temperature of the node
at the beginning of dt, i.e., TP > T 0S;p .
As shown in Fig. 4, at each time step dt, the procedure used to

Fig. 4. Flow chart of the computation procedure of auxiliary energy QA with on-off control.
 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906 897

simulate solar water heating systems with on-off control is  

computed in the following three steps: 4:69p rad 2p rad
d¼ sin ð284 þ nÞ :
36 365
1. Temperature TP is calculated using Eq. (5) with TI ¼ T 0S;b . If TI <
The solar radiation power flux reaching the collector may be
TP < T maxP ; entrance node p is evaluated according to the
approximated to the average solar radiation power during dt:
following criterion: p ¼ 1 if TP > T 0S;1 , or p ¼ kif k  2and
T 0S;k
1  TP > T 0S;k ; otherwise, i.e., if TP  TI or TP  T max _
P , CP ¼ _ G:
0 kW + C
1 . Gz
dt
2. If TB < TC , where TB ¼ T 0S;1 , flow rate C_ B ¼ C_ C , and auxiliary
Since the weather data is seldom available on time intervals of
energy QA is calculated using Eq. (6); otherwise, C_ B is calculated
less than one hour [2], the solar radiation power flux is normally
using Eq. (7), and QA ¼ 0 kWh.
approximated using dt ¼ 1 h.
3. Final temperatures T 1S;1 ; T 1S;2 ; …; T 1S;b are calculated using Eqs.
(8a)e(8e), since these will be initial temperatures
T 0S;1 ; T 0S;2 ; …; T 0S;b of the next time step. 2.6. Mains water temperature

The model developed by Araújo and Pereira [27] is used to

calculate mains temperature TM as a function of day n:
2.5. Solar radiation
 
2p rad
TM ¼ r1 T E þ r2 dTE expð
qÞsin n
q ;
The overall solar radiation energy (G, kWh m
2 ) reaching each 365
square meter of collector surface during dt is based on the HDKR
model [2,50,52], and its computational implementation is where T E (+ C) is the yearly mean outdoor temperature, and dTE (+ C)
described by Araújo and Pereira [27,28]. The following inputs are is the difference between maximum and minimum monthly out-
required: time step dt, the initial solar time (t 0 , h), the day of the door temperature means. Parameters r1 ¼ 0:87, r2 ¼ 2:3, and
year (n), the latitude (f, rad), the collector tilt (b, rad), the collector q ¼ 2:0 rad were obtained by data fitting. Since there were no
azimuth (g, rad), and the extraterrestrial (G0;z , kWh m
2 ), beam mains water temperature data for Portugal, temperature data from
(Gb;z , kWh m
2 ), and diffuse (Gd;z , kWh m
2 ) solar energy reaching six Spanish cities near the Portuguese border [53,54] were
each square meter of a horizontal surface during dt. The value of G is employed in the data fitting procedure.
computed as follows:
  3. Simulation parameters
G ¼ Gb;z þ Gd;z a Rb
1 þ cosb  b 3.1. Weather data
þGd;z ð1
aÞ 1 þ f sin3
2 2
  1
cosb Values of outdoor temperature TE and solar radiation fluxes G0;z ,
þ Gb;z þ Gd;z rg ; Gb;z , and Gd;z were obtained from files generated using computer
2
application CLIMAS-SCE, available for download from the web page
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
where a ¼ Gb;z =G0;z is the anisotropy index, f ¼ Gb;z =ðGb;z þ Gd;z Þ of the Portuguese national laboratory for energy and geology [55].
is the modulating factor, rg is the reflectivity of the surrounding These files are in the EnergyPlus weather format [56] and contain
surfaces, and Rb ¼ Sb =Sb;z is the ratio of beam radiation on the 8760 lines of hourly weather data for different Portuguese loca-
collector surface to that on a horizontal surface. Collector beam tions, i.e., day n ¼ 1; 2; …; 365, and initial time t 0 ¼ 0; 1; …; 23.
radiation is proportional to A single geographical location was chosen for all simulations:
municipality of Lisboa, Portugal, whose altitude is 109 m, longitude

Sb ¼ u1
u0 sindðsinfcosb
cosfsinbcosgÞ is
0:051p rad, and latitude f ¼ 0:216p rad.
 For the calculation of unconditioned temperature TU using Eq.
þ sinu1
sinu0 cosdðcosfcosb þ sinfsinbcosgÞ (2), a temperature TN ¼ 20 + C is a typical design value for indoor

conditions [48].

cosu1
cosu0 cosdsinbsing;

3.2. Design parameters

and horizontal beam radiation is proportional to
  Design parameters are properties or dimensions of the various
Sb;z ¼ u1
u0 sindsinf þ sinu1
sinu0 cosdcosf:
elements that form a solar water heating system, i.e., characteristics
of the solar collector, control scheme, storage tank, and working
However, if Rb < 0 (the sun is behind the collector), then Rb ¼ 0.
fluid. These parameters can be more or less easily modified
Angles u0 and u1 (rad) are the time angles at solar times t 0 and
accordingly in order to improve the energetic performance of the
t 1 ¼ t 0 þ dt, respectively. Time angle
solar heating system.
p rad As shown in Table 1, some design parameters were set as fixed
u¼ ðt
12 hÞ: values, while the others were allowed to vary in order to evaluate
12 h
their effect on the long-term performance of the solar heating
However, if u0 < usr < u1 (sunrise is within current time step), systems. With respect to the fixed parameters, maximum temper-
then u0 ¼ usr ; if u0 < uss < u1 (sunset is within current time step), ature TP;max is 5 + C below the temperature of boiling water at at-
then u1 ¼ uss . Angles usr and uss (rad) are the time angles at mospheric pressure; parameters h0 and FU are characteristic
sunrise and sunset, respectively: parameters of unglazed (UG) and single-glazed (SG) collectors [49];
collector tilt b equals latitude f, since it maximizes long-term
uss ¼
usr ¼ arctanð
tanftandÞ
absorbed solar radiation [2]; collector azimuth g corresponds to a
Angle d (rad) is the solar declination: south-oriented collector; reflectivity rg is a typical value for
898 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906

Table 1
Design parameters.

Element Parameter

Control parameters V_ P (production volume flow rate) 0.01e1 m3 h
1

TP;max (maximum production temperature) 95  C

Solar collector A (collector area) 0.2e20 m2

h0 (optical efficiency) 0.9 (UG), 0.8 (SG)
FU (heat loss coefficient) 0.02 (UG), 7  10
3 kW m
2 + C
1 (SG)
b (collector tilt) 0:216p rad
g (collector azimuth) 0 rad
rg (reflectivity of surrounding surfaces) 0.2

Storage tank VS (storage volume) 0.1e2 m3

US (storage heat transfer coefficient) 0 kW m
2 + C
1

Working fluid r (density) 988 kg m
3

c (specific heat capacity) 1:16  10
3 kWh kg
1  C
1

snowless surfaces [2]; coefficient US corresponds to a perfectly

thermally insulated storage tank (this value was used in all simu- rS DS
lS ¼ :
lations, except those in Section 4.5); fluid properties r and c are for b
water at atmospheric pressure and 50 + C [49].
It is assumed that the storage tank is a cylinder with a length-to-
diameter ratio rS ¼ 3, and all storage nodes k ¼ 1; 2; …; b are of
equal volume, i.e., VS;k ¼ VS =b. Therefore, depending on the num- 3.3. Consumption profile
ber of nodes b, the tank wall area of each node k can be calculated as
follows: The 2006 Portuguese regulation for the thermal performance of
buildings [57] assumes the following daily reference values of hot
water consumption: a volume vC ¼ 0:04 m3 per consumer at tem-
perature TC ¼ 60 + C.
8  
DS The so-called RAND profile [58], was used to model the daily
>
> pD þ l if b ¼ 1
>
> S
2 S distribution of hot water consumption. The RAND profile defines
>
<   the hourly fraction (xC ) of hot water that is consumed relative to the
AS;k ¼ DS ; total daily consumption.
>
> pDS þ lS if k ¼ 1; b and b  2
>
> 4 The consumption volume flow rate was approximated to the
>
:
pDS lS if 2  k  b
1 and b  2 average flow rate consumed during every hour:

x p v
where V_ C z C C C ;
1h

where pC is the number of consumers. For all simulations, it was

 1=3 assumed that pC ¼ 4, resulting in the volume flow rate profile
4VS shown in Fig. 5.
DS ¼ ;
prS

and 3.4. Performance output parameter

As stated in Section 2, yearly auxiliary energy SQA is the main

long-term performance simulation output from the solar water
heating models. However, as SQA is an extensive parameter that
depends on the yearly consumption of thermal energy (SQC , kWh),
the yearly solar fraction (F), an intensive and dimensionless
parameter defined as the ratio of the useful solar energy to the
thermal consumption load [2,50], was used to characterize the
performance of the simulated solar thermal systems:

SQC
SQA
F¼ :
SQC
The yearly consumption of thermal energy can be computed as
follows:

X
365
SQC ¼ CC ðTC
TM Þ;
n¼1

where CC ¼ rpC vC c, and pC vC ¼ 0:48 m3 represents the daily vol-

Fig. 5. Daily variation of consumption flow rate V_ C (based on the RAND profile [58]). ume of hot water consumption (Section 3.3).
 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906 899

Fig. 6. Variation of error dF with time step dt, where b ¼ 3; 4; 5, V_ P ¼ 0:05 m3 h
1 , Fig. 7. Variation of additional second-year days N
365 with convergence parameter
A ¼ 2; 4; 8 m2 , and VS ¼ 0:3; 0:4; 0:5 m3 . Legend: collector type. dT ref _ 3
1 2 3
S , where b ¼ 3; 4; 5, V P ¼ 0:05 m h , A ¼ 2; 4; 8 m , and VS ¼ 0:3; 0:4; 0:5 m .
Legend: collector type.

3.5. Numerical parameters

 
F ref
F 
Two important numerical parameters directly affect the accu-
dF ¼  
;
racy of the output results from the solar water heating simulations:  F ref 
time step dt and the convergence error of the storage temperatures
between two consecutive years. The first parameter is discussed in where F is the solar fraction at time step dt, F ref is the solar fraction
the next section, and the second in Section 3.5.2. at dt ref , and dt ref ¼ 10
4 h (10000 divisions per hour).
Another important numerical parameter is the number of nodes As shown in Fig. 6, for both UG and SG collector types, error dF
b assumed for the storage tank. It indicates the degree of stratifi- and the rate of increase of dF decrease with decreasing dt, so that dF
cation of the storage fluid between a fully-mixed tank, modeled by becomes practically zero when dt(0:01 h for all simulated condi-
a one-node tank, and a fully-stratified tank, modeled by a tank with tions. A value of dt ¼ 0:05 h was chosen for all simulations that
an infinite number of nodes [29]. The effect of this parameter is follow, since this is a good compromise between accuracy
evaluated in the sensitivity analysis presented in Section 4.2. (dF(0:15 %) and computational size (20 divisions per hour).
However, a time step of 0:1 h could also have been used without a
significant loss in accuracy (dF(0:3 %) and about the double of the
computation speed.

3.5.1. Time-step dependency 3.5.2. Stopping criterion for periodic convergence

As stated in Section 2.4, the convergence of the forward Euler At the beginning of each simulation (i.e., when n ¼ 1 and t 0 ¼
method, used to derive the explicit solutions (8a)e(8e) for storage 0 h), it is assumed that storage temperatures T 0S;1 ; T 0S;2 ; …; T 0S;b are
temperatures T 1S;1 ;T 1S;2 ;…;T 1S;b , is dependent on the size of time step all equal to mains water temperature TM . However, one year ahead
dt, being restricted to short time steps. It is also known that, (i.e., when n ¼ 366 and t 0 ¼ 0 h), temperatures T 0S;1 ; T 0S;2 ; …; T 0S;b
generally, the error committed by the Euler approximation is pro- will be different from TM . Therefore, in order to obtain a yearly
portional to the time step [51]. Therefore, an optimum value of dt periodic simulation, the simulation has to continue during the
must be small enough in order to guarantee convergence and suf- beginning of the second simulation year until the following crite-
ficiently accurate solutions, but not so small that the computation rion is satisfied:
process becomes prohibitively slow.
Solar thermal simulations were performed for the following 
max dTS;1 ; dTS;2 ; …; dTS;b  dT ref
S :
time steps: dt ¼ 2:5  10
4 ; 5  10
4 ; 10
3 ; 2:5  10
3 ; …; 1 h
(from 4000 to one divisions per hour). The change in time step only Series dTS;1 ; dTS;2 ; …; dTS;b is formed by the relative difference
affects the numerical model developed for the thermal system between the storage temperature on day n (T new S;k ) and the tem-
(Section 2.4). For the radiation model presented in Section 2.5, the perature on day n
365 (T old 0
S;k ) at the same time t and nodes k ¼
weather data (Section 3.1), and the consumption profile (Section 1; 2; …; b:
3.3), time step dt ¼ 1 h. The simulated design parameters were  new 
T old 
chosen so as not to be far from the range of the typical parameters  S;k
T S;k 
dTS;k ¼  :
used with the consumption load defined in Section 3.3: number of  T new
S;k

storage nodes b ¼ 3; 4; 5, production flow rate V_ P ¼ 0:05 m3 h
1 ,
collector area A ¼ 2; 4; 8 m2 , and storage volume VS ¼ 0:3; 0:4; The value of dT ref
S controls the convergence of the storage
0:5 m3 . temperatures on the second simulation year with respect to the
The relative error (dF) of the solar fraction at a given time step dt temperatures on the first year.
from the solar fraction at a reference time step (dt ref ) was used to Fig. 7 shows the variation of the number of additional second-
evaluate the convergence of the simulated outputs with respect to year days (N
365, where N is the total number of simulation
decreasing time steps: days) with convergence parameter dT ref S . The same design
900 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906


Fig. 9. Variation of F max and V_ P with number of nodes b for the UG collector. Legend:

Fig. 8. Variation of solar fraction F with flow rate V_ P , where b ¼ 3; 4; 5, A ¼ 2; 4; 8 m2 , A, VS . (a) Maximum solar fraction F max . (b) Optimum flow rate V_ P .
and VS ¼ 0:3; 0:4; 0:5 m3 . Legend: A. (a) Unglazed collector. (b) Single-glazed collector.

Therefore, for different design parameters, an optimum value of


parameters used for the time-step dependency study (Section production flow rate (V_ P , m3 h
1 ) should be searched so as to
3.5.1) were selected here for both UG and SG collector types, i.e., maximize F, i.e.,
b ¼ 3; 4; 5, V_ P ¼ 0:05 m3 h
1 , A ¼ 2; 4; 8 m2 , and VS ¼ 0:3; 0:4;
0:5 m3 . For the stopping criterion of all simulations, the following F max ¼ maxF;
V_ P
convergence parameter was selected: dT ref S ¼ 10
3 , since N
365(
15 (less than 9600 additional time steps for dt ¼ 0:05 h), a good
and
compromise between accuracy and computation speed.

V_ P ¼ argmaxF;
V_P
4. Sensitivity analysis
max
subjected to 0 m3 h
1  V_ P  V_ P . Even though Araújo and Per-
4.1. Production flow rate eira [28] limited the maximum flow rate of proportional control
systems to 1 m3 h
1 , for the present system, F was found to be
According to Eq. (5), production temperature TP is strongly almost unaffected by the variation of V_ P above 0:5 m3 h
1 , so that
dependent on flow rate C_P , and hence on V_P . Therefore, as the value max
V_ P ¼ 0:5 m3 h
1 for all subsequent simulations. It is worth
max
of TP regulates whether the production circuit is either turned on or noting that the value of V_ P must be high enough to guarantee
off, the preset value of V_P has a strong impact on the long-term that F max max
is within flow rates of 0eV_ P , but not too high, since the
max
performance of the solar systems. time to compute F max increases with increasing V_ P .
Fig. 8 shows the relationship between yearly solar fraction F and Package Optim, a package for univariate and multivariate opti-
V_ P (0.01e0:5 m3 h
1 ) for the UG and SG collector types and the mization in Julia [59,60], was used to find maximum solar fraction

following conditions: number of nodes b ¼ 3; 4; 5, collector area F max and optimum flow rate V_ P by means of Brent’s method [61].
A ¼ 2; 4; 8 m2 , and storage volume VS ¼ 0:3; 0:4; 0:5 m3 . Solar
fraction F increases very rapidly with increasing V_ P for low values of
V_ P , reaching a maximum value (F max ) from which F decreases or 4.2. Degree of stratification
remains approximately constant with further increasing V_ P . The
value of F max and the value of V_ P for which F max is reached is The number of nodes b regulates the degree of stratification
strongly dependent on the size of collector area A, but it also de- achievable by the storage tank: the fluid in a storage tank with one
pends on the other design parameters. node does not stratify (i.e., it is fully mixed), the degree of
 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906 901


 Fig. 11. Variation of F max and V_ P with collector area A for the UG collector. Legend: b,
Fig. 10. Variation of F max and V_ P with number of nodes b for the SG collector. Legend: 
 VS . (a) Maximum solar fraction F max . (b) Optimum flow rate V_ P .
A, VS . (a) Maximum solar fraction F max . (b) Optimum flow rate V_ P .

stratification increases as b increases, and, as b tends to infinity, the keeps decreasing with increasing b. The reason for the low values of

tank becomes fully stratified [2,29]. V_ P when b is high is that when thermal stratification becomes
Figs. 9a and 10a show the variation of maximum solar fraction significant, the fluid must enter the storage tank at low speeds in
F max with b (1  b  10) for the UG and SG collector types, order to avoid the destruction of stratification [2,31,39].
respectively, and the following design conditions: collector area Duffie and Beckman [2] suggested that three to four nodes
A ¼ 2; 4; 8 m2 , and storage volume VS ¼ 0:3; 0:4 m3 . For the UG represents a reasonable compromise between an unstratified sys-
collector, F max increases 25e28 % when b increases from 1 (fully tem and the limiting situation of full stratification, and, therefore, a
mixed) to 4 (a realistic degree of stratification [2]), which is a sig- value of b ¼ 4 was chosen as the stratification level for all subse-
nificant improvement in system performance. For the SG collector, quent simulations.
F max increases a little less (5e16 %) when b increases from 1 to 4.
The increase in F max is higher for the lowest collector area (i.e., 4.3. Collector area
when A ¼ 2 m2 ). The reason for the higher performance
improvement achieved by the UG collector is due to the lower value In addition to efficiency parameters h0 and FU found in Eq. (4),
of F max , since, as the solar fraction is limited to unity, there is more collector area A is the main design parameter, as the performance of
room for solar fraction improvement when the solar fraction is low. solar thermal systems is much more sensitive to A than to any other
The increase in the simulated values of F max due to thermal strat- parameter [2].
ification is in line with the findings reported in the literature Figs. 11a and 12a show the variation of maximum solar fraction
(Section 1): system efficiency rises in the range of 5e20 % have F max with A (0.2e20 m2 ) for the UG and SG collectors, respectively,
been reported in the literature [34e40], although efficiency im- and the following conditions: number of nodes b ¼ 1; 4, and stor-
provements as high as 38 % have also been reported by Hollands age volume VS ¼ 0:3; 0:4; 0:5 m3 . As expected, for both collector
and Lightstone [41]. types, solar fraction F max increases with increasing A, but its rate of

The values of optimum flow rate V_ P are plotted in Figs. 9b and increase decreases with increasing A, especially when F max is closer
10b against b. When the number of nodes is low (b ¼ 1; 2 for the to unity. Solar fraction F max is significantly higher when b ¼ 4,

UG collector, and b ¼ 1; 2; 3 for the SG collector), flow rate V_ P ¼ confirming the improvement of system performance due to the
max
_
VP 3
1
¼ 0:5 m h , indicating that solar fraction F increases stratification of the storage fluid, as discussed in Section 4.2.

monotonically with V_ P in the range 0e0:5 m3 h
1 . However, for Figs. 11b and 12b show the variation of optimum flow rate V_ P
higher values of b (b  3 for the UG collector, and b  4 for the SG with A for the UG and SG collectors, respectively, when b ¼ 4. For

 
collector), flow rate V_ P drops down to values of V_ P (0:1 m3 h
1 and the UG collector, flow rate V_ P increases approximately from 0.02 to

1
0:06 m h as A increases from 0.2 to 20 m2 , but its rate of increase
3
902 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906


Fig. 12. Variation of F max and V_ P with collector area A for the SG collector. Legend: b,

VS . (a) Maximum solar fraction F max . (b) Optimum flow rate V_ P .
Fig. 13. Variation of solar fraction F with collector area A and production flow rate V_ P .
Legend: b, V_ P . (a) Unglazed collector. (b) Single-glazed collector.

decreases with increasing A. Likewise, for the SG collector, V_ P z

3
1 _
0:02 m h for the lowest values of A, V P increases with increasing
 and the value of A for which F is maximum both increase with
A, but, for Aa4 m2, V_ P increases very quickly with A, reaching flow increasing V_ P . This behavior is justified by the relationship between
rates of approximately 0:5 m3 h
1 for the highest values of A. temperature TP and dependent variables A, h0 , FU, and V_ P , as

Although not shown in Figs. 11b and 12b, when b ¼ 1, V_ P z established in Eq. (5) and discussed in the previous paragraph,
0:5 m3 h
1 for the whole range of simulated values of A. which also clarifies the reason for the high values of optimum flow
Considering collector areas in the range of 2e8 m2 , optimum 
rate V_ P shown in Figs. 11b and 12b when A is high. The unstratified
values of production flow rate per square meter of collector area (b ¼ 1) values presented in Fig. 13 agree with the values presented

(V_ P =A) were found to be in the range of 0.006e0:016 m h
1 and by Araújo and Pereira [27] for a similar collector type and similar
0.01e0:014 m h
1 for the UG and SG collectors, respectively, where design conditions.

V_ P =A decreases with increasing A. In Particular, when A ¼ 4 m2 and The improvement of the energy performance of stratified sys-
 
VS ¼ 0:4 m3 ; V_ P =A ¼ 0:01 m h
1 for the UG collector, and V_ P = A ¼ tems may be also appreciated by the fact that the value of

0:011 m h for the SG collector. These simulated values of V_ P = A

1
maximum F and the corresponding value of V_ P both increase
are in close agreement with those found in the literature (i.e., in the significantly when b increases from 1 to 4, as shown in Fig. 13.
order of 0:01 m h
1 ) [2,41,42].
According to Eq. (5), production temperature TP increases line-
arly with increasing A. Therefore, the production circuit is generally 4.4. Storage volume
turned off more frequently and for longer periods when A is high, as
TP exceeds maximum temperature T max P more often, especially The variation of maximum solar fraction F max with storage
when collector efficiency parameters h0 and FU and solar radiation volume VS (0.1e2 m3 ) is shown in Figs. 14a and 15a for the UG and
G_ are also high, contributing to a reduction in yearly solar fraction F. SG collectors, respectively, and the following conditions: number of
Conversely, according to Eq. (5), increasing production flow rate C_ P , nodes b ¼ 1; 4, and collector area A ¼ 2; 4; 8 m2 . The effect of
and hence V_P ; has the effect of decreasing temperature TP . storage volume on solar fraction is not very significant: solar frac-
Fig. 13 shows the variation of solar fraction F with A for the UG tion F max increases slightly with increasing volume VS when VS is
and SG collectors and fixed values of flow rate V_ P , where b ¼ 1; 4, low, tending to an almost constant value for higher values of VS .
V_ P ¼ 0:02;0:05;0:1 m3 h
1 , and VS ¼ 0:4 m3 . For low values of V_ P , F As shown in Figs. 14b and 15b, when b ¼ 4, optimum flow rate

increases with increasing A when A is low, but, as A increases, V_ P (0:04 m3 h
1 and remains more or less unchanged with vary-
especially for the most efficient SG collector, F reaches a maximum ing VS for the two lowest collector areas (i.e., A ¼ 2; 4 m2 ); when

and starts decreasing with increasing A. The value of maximum F A ¼ 8 m2 , V_ P increases with increasing VS up to approximately
 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906 903

 
Fig. 14. Variation of F max and V_ P with storage volume VS for the UG collector. Legend: Fig. 15. Variation of F max and V_ P with storage volume VS for the SG collector. Legend:
 
b, A. (a) Maximum solar fraction F max . (b) Optimum flow rate V_ P . b, A. (a) Maximum solar fraction F max . (b) Optimum flow rate V_ P .


0:06 m3 h
1 for the UG collector, but V_ P has a somewhat erratic insulation layer is considered), respectively, and hS (kW m
2 + C
1 )
variation with varying VS for the SG collector. Although not shown is the coefficient of convection on the external surfaces of the tank.

in Figs. 14b and 15b, when b ¼ 1, V_ P z0:5 m3 h
1 , the maximum However, due to thermal leaks through tank connections, this
allowed value of V_ P , for all simulated values of VS . equation normally underestimates the value of US [2], and, there-
fore, slightly overestimated values were assumed for the thermal
4.5. Storage tank losses conductivity and coefficient of convection [62]: kS ¼ 5 
10
5 kW m
1 + C
1 , and hS ¼ 0:01 kW m
2 + C
1 .
In Sections 4.1e4.4, all simulations were performed assuming a Similarly to Figs. 9a, 10a and 16 shows the variation of maximum
perfectly thermally insulated storage tank. Araújo and Pereira solar fraction F max with b (1  b  10) for the UG and SG collector
[27,28] concluded that tank thermal losses can be ignored above types and the following conditions: collector area A ¼ 2; 4; 8 m2 ,
insulation thicknesses of the order of 0:2 m. However, in practice, and storage volume VS ¼ 0:4 m3 . However, in addition to the
storage tanks are often covered with thinner layers of insulation. In perfectly insulated tank (US ¼ 0 kW m
2 + C
1 ), Fig. 16 also shows
addition, the simulated solar fractions presented in Sections sec: the values of F max for tank wall thickness dxS ¼ 0:05 m (US ¼ 9:1 
flow rate,sec: strat,sec: coll area,sec: stor vol seem to be excessively 10
4 kW m
2 + C
1 ), indicating a considerable reduction in F max
high regarding the solar fractions normally reported in the litera- for the less efficient tank. For example, with the SG collector, when
ture for comparable systems. For example, with the SG collector, VS ¼ 0:4 m3 and A ¼ 8 m2 , F max ¼ 0:88and 0.93 for b ¼ 1 and 4,
when VS ¼ 0:4 m3 and A ¼ 8 m2 , F max ¼ 0:92 and 0.97 for b ¼ 1 respectively; when A ¼ 4 m2 , F max ¼ 0:75and 0.85 for b ¼ 1 and 4,
and 4, respectively; when A ¼ 4 m2 , F max ¼ 0:83 and 0.91 for b ¼ 1 respectively. These values are in line with those found in similar
and 4, respectively. As a result, this section investigates the effect of solar water heating systems, since, considering the collector type,
tank thermal losses on the overall performance of solar water the consumption load, and the subtropical climate of Lisboa (about
heating systems. 1800 kWh m
2 of global horizontal solar irradiation per year [63]),
The heat transfer coefficient of the tank walls can be estimated an area of 4 m2 should suffice to provide a solar fraction of the order
as follows [62]: of 0.8 [64].
For wall thickness dxS ¼ 0:05 m, it was found that F max in-
1 creases 25e28 % and 6e20 % for the UG and SG collectors, respec-
US ¼ ; tively, when b increases from 1 to 4. Even though there is a small
dxS =kS þ 1=hS
improvement in the increase of F max for the SG collector, these
where dxS (m) and kS (kW m
1 + C
1 ) are the effective thickness values are very similar to those reported in Section 4.2 for the
and thermal conductivity of the tank walls (generally, only the perfectly insulated tank.
904 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906

Fig. 16. Variation of maximum solar fraction F max with number of nodes b for a Fig. 17. Variation of relative maximum solar fraction difference dF max with insulation
perfectly insulated tank (US ¼ 0 kW m
2 + C
1 ) and a tank with wall thickness dxS ¼ thickness dxS , where VS ¼ 0:4 m3 . Legend: b, A. (a) Unglazed collector. (b) Single-
0:05 m (US ¼ 9:1  10
4 kW m
2 + C
1 ), where VS ¼ 0:4 m3 . Legend: A, US . (a) Un- glazed collector.
glazed collector. (b) Single-glazed collector.

5. Conclusions
The relative decrease in maximum solar fraction (dF max ) was
used to quantify the effect of different insulation thicknesses on
For the estimation of the long-term performance of solar water
system performance:
heating systems, there is generally no reason to produce overly
complex models, as these tend to be computationally inefficient,
and the apparent accuracy added by fine system-detail modeling is
mostly dissipated by the errors introduced with the assumed
F max;0
F max weather data and consumption profile. The present work presents a
dF max ¼ ; simplified solar thermal model, which promotes linearity and ex-
F max;0
cludes implicit solutions, for the rapid computation of the long-
where F max;0 is the maximum solar fraction for a perfectly insulated term performance, i.e., yearly solar fraction F, of direct solar water
tank, and F max is the maximum solar fraction for a given level of heating systems using the on-off control scheme. Thermal strati-
storage losses. fication was also included by means of a simple one-dimensional
Fig. 17 shows the variation of relative difference dF max with multi-node model, in which the storage tank is divided into b
thickness dxS (0e0:3 m) for the UG and SG collector types and the vertically stacked nodes. Number of nodes b indicates the degree of
following conditions: number of nodes b ¼ 1; 4, collector area A ¼ stratification achievable by the storage tank.
2; 4; 8 m2 , and storage volume VS ¼ 0:4 m3 . The value of dF max and Model simulations were performed for the evaluation of the
its rate of decrease both decrease with increasing thickness dxS , so relative impact of the most significant system parameters (i.e., time
that dF max tends to zero as dxS becomes very large. When dxS ¼ step dt, number of nodes b, production flow rate V_ P , solar collector
0:2 m (US ¼ 2:4  10
4 kW m
2 + C
1 ), dF max is in the range of area A, storage tank volume VS , and storage heat transfer coefficient
2.1e3:4 % for the UG collector, and dF max is in the range of 1e5:2 % US ) on the yearly solar fraction, enabling the model to be evaluated
for the SG collector. The value of dF max is generally higher for the using published results from comparable systems. The Portuguese
fully mixed case (b ¼ 1), and it increases with decreasing A. For municipality of Lisboa, with a subtropical climate, was used for the
lower values of dxS , the value of dF max becomes too high for the weather database, two different solar collector types were selected,
storage losses to be neglected (e.g., when dxS ¼ 0:05 m, dF max is in and the walls of the storage tank were assumed to be perfectly
the range of 7.3e11:2 % for the UG collector, and it is in the range of thermally insulated.
3.5e16:4 % for the SG collector). A time-step dependency analysis was performed by evaluating
 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906 905

the convergence of solar fraction F with decreasing time step dt. A References
good compromise between accuracy (a maximum error of about
[1] F. Mauthner, W. Weiss, M. Spo € rk-Dür, Solar Heat Worldwide: Markets and
0:15 %) and computation speed (20 divisions per hour) was found
Contribution to the Energy Supply 2014, IEA Solar Heating & Cooling Pro-
for dt ¼ 0:05 h, and, thus, this time step was used for all subsequent gramme, Gleisdorf, Austria, 2016.
simulations. [2] J.A. Duffie, W.A. Beckman, Solar Engineering of Thermal Processes, fourth ed.,
The dynamics of solar water heating systems is one-year peri- Wiley, Hoboken, NJ, 2013 https://doi.org/10.1002/9781118671603.
[3] Z. Wang, W. Yang, F. Qiu, X. Zhang, X. Zhao, Solar water heating: from theory,
odic. In order to achieve periodic convergence in the simulations, a application, marketing and research, Renew. Sustain. Energy Rev. 41 (2015)
stopping criterion was chosen such that the relative difference of 68e84, https://doi.org/10.1016/j.rser.2014.08.026.
the storage temperatures between two consecutive simulation [4] S.R. Schiller, M.L. Warren, D.M. Auslander, Comparison of proportional and on/
off solar collector loop control strategies using a dynamic collector model,
years should not exceed 10
3 , resulting in about 15 additional J. Sol. Energy Eng. 102 (4) (1980) 257e262, https://doi.org/10.1115/
simulation days after the first simulation year. 1.3266189.
It was found that solar fraction F increases with increasing [5] M. Boji
c, S. Kalogirou, K. Petronijevi
c, Simulation of a solar domestic water
production flow rate V_ P when V_ P is low. For fully-mixed storage, heating system using a time marching model, Renew. Energy 27 (3) (2002)
441e452, https://doi.org/10.1016/S0960-1481(01)00098-2.
the rate of increase of F decreases with increasing V_ P , tending to an [6] A. Georgiev, Simulation and experimental results of a vacuum solar collector
almost constant value when V_ P is high. However, when tank system with storage, Energy Convers. Manag. 46 (9e10) (2005) 1423e1442,
https://doi.org/10.1016/j.enconman.2004.07.002.
stratification is included, particularly when number of nodes ba3,
[7] R. Kumar, M.A. Rosen, Thermal performance of integrated collector storage
the solar fraction reaches a maximum F max at an optimum flow rate solar water heater with corrugated absorber surface, Appl. Therm. Eng. 30 (13)

V_ P , from which F starts decreasing with further increasing V_ P . (2010) 1764e1768, https://doi.org/10.1016/j.applthermaleng.2010.04.007.
Therefore, for all subsequent simulations, a numerical procedure [8] Y.-D. Kim, K. Thu, H.K. Bhatia, C.S. Bhatia, K.C. Ng, Thermal analysis and per-
formance optimization of a solar hot water plant with economic evaluation,
was applied for finding maximum solar fraction F max and corre- Sol. Energy 86 (5) (2012) 1378e1395, https://doi.org/10.1016/

sponding optimum flow rate V_ P . j.solener.2012.01.030.
Maximum solar fraction F max increases with number of nodes b. [9] R. Kicsiny, Z. Varga, Real-time state observer design for solar thermal heating
systems, Appl. Math. Comput. 218 (23) (2012) 11558e11569, https://doi.org/
Increments in the range of 5e28 % were observed when b increases 10.1016/j.amc.2012.05.040.
from 1 to 4. These findings are in close agreement with results [10] R. Kicsiny, I. Farkas, Improved differential control for solar heating systems,
reported in the literature for similar systems. The increase in F max is Sol. Energy 86 (11) (2012) 3489e3498, https://doi.org/10.1016/
j.solener.2012.08.003.
superior for less efficient collectors and lower collector areas. [11] MathWorks-Makers of MATLAB and Simulink, 2019. https://www.
In low-stratified systems, especially when number of nodes mathworks.com/.

b(3, optimum flow rate V_ P is the maximum flow rate allowed by [12] TRNSYS-Official Website, 2019. http://sel.me.wisc.edu/trnsys/.
[13] J. Buz

, R. Ne
meth, Modelling and simulation aspects of a
the system (in the present work, 0:5 m3 h
1 ). However, in stratified as, I. Farkas, A. Biro
solar hot water system, Math. Comput. Simulat. 48 (1) (1998) 33e46, https://
systems, in order to avoid the destruction of stratification, optimum doi.org/10.1016/S0378-4754(98)00153-0.

flow rate V_ P is reduced to values of the order of 0:01 m3 h
1 per [14] V. Badescu, Optimal control of flow in solar collector systems with fully mixed
water storage tanks, Energy Convers. Manag. 49 (2) (2008) 169e184, https://
square meter of collector area. For example, when b ¼ 4 and the
 doi.org/10.1016/j.enconman.2007.06.022.
collector area is between 2 and 8 m2 , V_ P is in the range of [15] A. Hobbi, K. Siddiqui, Optimal design of a forced circulation solar water

0.006e0:016 m h per square meter of collector area, where V_ P is
3
1
heating system for a residential unit in cold climate using TRNSYS, Sol. Energy
higher for lower collector areas. These values are in line with results 83 (5) (2009) 700e714, https://doi.org/10.1016/j.solener.2008.10.018.
[16] L.M. Ayompe, A. Duffy, S.J. McCormack, M. Conlon, Validated TRNSYS model
reported in the literature for comparable systems. for forced circulation solar water heating systems with flat plate and heat pipe
Maximum solar fraction F max increases with increasing collector evacuated tube collectors, Appl. Therm. Eng. 31 (8e9) (2011) 1536e1542,
area A, but its rate of increase decreases with increasing A, so that https://doi.org/10.1016/j.applthermaleng.2011.01.046.
[17] M.C. Rodríguez-Hidalgo, P.A. Rodríguez-Aumente, A. Lecuona, M. Legrand,
F max tends to an almost constant value for high values of A. R. Ventas, Domestic hot water consumption vs. solar thermal energy storage:
The effect of storage volume VS on maximum solar fraction F max the optimum size of the storage tank, Appl. Energy 97 (2012) 897e906,
is not very significant, since F max increases only slightly with https://doi.org/10.1016/j.apenergy.2011.12.088.
[18] F. Mosallat, T. ElMekkawy, D.L. Friesen, T. Molinski, S. Loney, E.L. Bibeau,
increasing VS when VS is low, tending to an almost constant value Modeling, simulation and control of flat panel solar collectors with thermal
for higher values of VS . storage for heating and cooling applications, Procedia Comput. Sci. 19 (2013)
The effect of tank thermal losses on maximum solar fraction 686e693, https://doi.org/10.1016/j.procs.2013.06.091.
[19] M.J.R. Abdunnabia, K.M.A. Alakder, N.A. Alkishriwi, S.M. Abughres, Experi-
F max were analyzed through the variation of storage heat transfer
mental validation of forced circulation of solar water heating systems in
coefficient US . It was found that, unless the tank walls are covered TRNSYS, Energy Procedia 57 (2014) 2477e2486, https://doi.org/10.1016/
by a rather thick layer of thermal insulation (about 0:2 m), storage j.egypro.2014.10.257.
[20] J. Buza
s, R. Kicsiny, Transfer functions of solar collectors for dynamical anal-
tank losses cannot be ignored. For thinner layers of thermal insu-
ysis and control design, Renew. Energy 68 (2014) 146e155, https://doi.org/
lation (about 0:05 m), the simulated solar fractions are in good 10.1016/j.renene.2014.01.037.
agreement with published results. [21] R. Kicsiny, New delay differential equation models for heating systems with
pipes, Int. J. Heat Mass Transf. 79 (2014) 807e815, https://doi.org/10.1016/
j.ijheatmasstransfer.2014.08.058.
[22] R. Kicsiny, J. Nagy, C. Szalo
ki, Extended ordinary differential equation models
Author contributions section for solar heating systems with pipes, Appl. Energy 129 (2014) 166e176,
https://doi.org/10.1016/j.apenergy.2014.04.108.
[23] R. Kicsiny, Multiple linear regression based model for solar collectors, Sol.
Anto
nio Araújo: Conceptualization; Methodology; Software;
Energy 110 (2014) 496e506, https://doi.org/10.1016/j.solener.2014.10.003.
Validation; Formal analysis; Investigation; Writing e Original [24] R. Kicsiny, Transfer functions of solar heating systems for dynamic analysis
Draft; Visualization; Project administration. Rui Silva: Software; and control design, Renew. Energy 77 (2015) 64e78, https://doi.org/10.1016/
j.renene.2014.12.001.
Resources; Writing e Review & Editing. [25] C.N. Antoniadis, G. Martinopoulos, Optimization of a building integrated solar
thermal system with seasonal storage using TRNSYS, Renew. Energy 137
(2019) 56e66, https://doi.org/10.1016/j.renene.2018.03.074.
[26] K. Kubin
ski, Ł. Szabłowski, Dynamic model of solar heating plant with sea-
Declaration of competing interest
sonal thermal energy storage, Renew. Energy 145 (2020) 2025e2033, https://
doi.org/10.1016/j.renene.2019.07.120.
The authors declare that they have no known competing [27] A. Araújo, V. Pereira, Solar thermal modeling for rapid estimation of auxiliary
energy requirements in domestic hot water production: on-off flow rate
financial interests or personal relationships that could have
control, Energy 119 (2017) 637e651, https://doi.org/10.1016/
appeared to influence the work reported in this paper.
906 A. Araújo, R. Silva / Renewable Energy 150 (2020) 891e906

j.energy.2016.11.025. [44] S. Fertahi, A. Jamil, A. Benbassou, Review on solar thermal stratified storage
[28] A. Araújo, V. Pereira, Solar thermal modeling for rapid estimation of auxiliary tanks (STSST): insight on stratification studies and efficiency indicators, Sol.
energy requirements in domestic hot water production: proportional flow Energy 176 (2018) 126e145, https://doi.org/10.1016/j.solener.2018.10.028.
rate control, Energy 138 (2017) 668e681, https://doi.org/10.1016/ [45] Y.H. Zurigat, K.J. Maloney, A.J. Ghajar, A comparison study of one-dimensional
j.energy.2017.07.109. models for stratified thermal storage tanks, J. Sol. Energy Eng. Trans. ASME
[29] H.P. Garg, S.C. Mullick, A.K. Bhargave, Solar Thermal Energy Storage, D. Reidel, 111 (3) (1989) 204e210, https://doi.org/10.1115/1.3268308.
Dordrecht, Holland, 1985. [46] J. Bezanson, A. Edelman, S. Karpinski, V.B. Shah, Julia, A fresh approach to
[30] Y.P. Chandra, T. Matuska, Stratification analysis of domestic hot water storage numerical computing 59 (1) (2014) 65e98, https://doi.org/10.1137/
tanks: a comprehensive review, Energy Build. 187 (2019) 110e131, https:// 141000671.
doi.org/10.1016/j.enbuild.2019.01.052. [47] The Julia Language, 2019. https://julialang.org/.
[31] Y.M. Han, R.Z. Wang, Y.J. Dai, Thermal stratification within the water tank, [48] Residential cooling and heating load calculations, in: ASHRAE Hand-
Renew. Sustain. Energy Rev. 13 (5) (2009) 1014e1026, https://doi.org/ bookdFundamentals (SI Edition), American Society of Heating, Refrigerating,
10.1016/j.rser.2008.03.001. and Air-Conditioning Engineers, Atlanta, GA, 2009. Ch. 17, pp. 17.1e17.16.
[32] J. Fan, S. Furbo, Thermal stratification in a hot water tank established by heat [49] B. Bourges, European Simplified Methods for Active Solar System Design,
loss from the tank, Sol. Energy 86 (11) (2012) 3460e3469, https://doi.org/ Kluwer Academic Publishers, Dordrecht, Holland, 1991.
10.1016/j.solener.2012.07.026. [50] S.A. Kalogirou, Solar Energy Engineering: Processes and Systems, second ed.,
[33] A. Li, F. Cao, W. Zhang, B. Shi, H. Li, Effects of different thermal storage tank Academic Press, Oxford, England, 2014.
structures on temperature stratification and thermal efficiency during [51] K.E. Atkinson, W. Han, D. Stewart, Numerical Solution of Ordinary Differential
charging, Sol. Energy 173 (2018) 882e892, https://doi.org/10.1016/ Equations, Wiley, Hoboken, NJ, 2011, https://doi.org/10.1002/
j.solener.2018.08.025. 9781118164495.
[34] M.K. Sharp, R.I. Loehrket, Stratified thermal storage in residential solar energy [52] D.T. Reindl, W.A. Beckman, J.A. Duffie, Evaluation of hourly tilted surface ra-
applications, J. Energy 3 (2) (1979) 106e113. diation models, Sol. Energy 45 (1) (1990) 9e17, https://doi.org/10.1016/0038-
[35] C.W.J. van Koppen, J.P.S. Thomas, W.B. Veltkamp, The actual benefits of 092X(90)90061-G.
thermally stratified storage in a small and medium size solar systems, in: [53] UNE 94002, Instalaciones Solares Te
rmicas para Produccio
n de Agua Caliente
Proceedings of the ISES Solar World Congress, Atlanta, GA, 1979, pp. 576e580. SanitariadC
alculo de la Demanda de Energía Te
rmica, Asociacio
n Espan
~ ola de
[36] W.B. Veltkamp, Thermal stratification in heat storages, in: C. den Ouden (Ed.), Normalizacio
n y Certificacio
n, Madrid, Spain, 2005.
Thermal Storage of Solar Energy, Martinus Nijhoff Publishers, Dordrecht, [54] Guía Te
cnicadCondiciones Clim
aticas Exteriores de Proyecto, Instituto para la
Holland, 1981, pp. 47e59. Diversificacio
n y Ahorro de la Energía, Madrid, Spain, 2010.
[37] R.L. Cole, F.O. Bellinger, Thermally stratified tanks, ASHRAE Transact. 88 (2) [55] LNEGdLaborato
rio Nacional de Energia e Geologia, 2019. http://www.lneg.pt/
(1982) 1005e1017. .
[38] M.D. Wuestling, S.A. Klein, J.A. Duffie, Promising control alternatives for solar [56] EnergyPlus, 2019. https://energyplus.net/.
water heating systems, J. Sol. Energy Eng. Trans. ASME 107 (3) (1985) [57] S. Camelo, C.P. Santos, A.
Ramalho, C. Horta, H. Gonçalves, E. Maldonado,
215e221, https://doi.org/10.1115/1.3267681. Manual de Apoio a  Aplicaç~ao do RCCTE, INETI, Lisboa, Portugal, 2006.
[39] J.A. Duffie, Overview of modelling and simulation, in: G. Lo € f (Ed.), Active Solar [58] J.J. Mutch, Residential Water Heating: Fuel Conservation, Economics, and
Systems, MIT Press, Cambridge, MA, 1993, pp. 3e18. Ch. 1. Public Policy (R-1498-NSF), RAND Corporation, Santa Monica, CA, 1974.
[40] N.K. Ghaddar, Stratified storage tank influence on performance of solar water [59] P.K. Mogensen, A.N. Riseth, Optim: a mathematical optimization package for
heating system tested in Beirut, Renew. Energy 4 (8) (1994) 911e925, https:// Julia, J. Open Source Softw. 3 (24) (2018) 615, https://doi.org/10.21105/
doi.org/10.1016/1352-0237(94)90033-7. joss.00615.
[41] K.G. Hollands, M.F. Lightstone, A review of low-flow, stratified-tank solar [60] JuliaNLSolvers-GitHub, 2019. https://github.com/JuliaNLSolvers.
water heating systems, Sol. Energy 43 (2) (1989) 97e105, https://doi.org/ [61] R.P. Brent, Algorithms for Minimization without Derivatives, Prentice-Hall,
10.1016/0038-092X(89)90151-5. Englewood Cliffs, NJ, 1973.
[42] M. Spanggaard, T. Schwaner, Basics of Thermal Stratification in Solar Thermal [62] F.P. Incropera, D.P. DeWitt, T.L. Bergman, A.S. Lavine, Fundamentals of Heat
Systems, Tech. Rep., EyeCular Technologies, Copenhagen, Denmark, 2015. and Mass Transfer, sixth ed., Wiley, Hoboken, NJ, 2006.
[43] M.Y. Haller, C.A. Cruickshank, W. Streicher, S.J. Harrison, E. Andersen, S. Furbo, [63] R. Vieira, J. Neto, A.
Silva, Radiaç~ao Solar Global AnualdPortugal Continental
Methods to determine stratification efficiency of thermal energy storage 2002e2007, Tech. Rep, Instituto de Meteorologia, Lisboa, Portugal, 2008.
processesdreview and theoretical comparison, Sol. Energy 83 (10) (2009) [64] Planning and Installing Solar Thermal Systems: A Guide for Installers, Archi-
1847e1860, https://doi.org/10.1016/j.solener.2009.06.019. tects and Engineers, second ed., Earthscan, London, England, 2010.